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ABSTRACT 


STRUCTURAL  OPTIMIZATION  INCLUDING 
CENTRIFUGAL  AND  CORIOLIS  EFFECTS 


Howard  Dwight  Gans 


Chairman:  William  J.  Anderson 


-'“  This  dissertation,  investigates  the  effects  of  centrifugal  and  Coriolis  forces  on  the 
mode  shapes  and  frequencies  of  a  rotating  system/^Optimal  redesign  is  then  done  to 
favorably  alter  the  mod The  rotational  effects  have  a  profound  influence  on  the 
eigenfrequencies;  this  is  important  in  optimal  structural  redesign  where  the  frequencies 
must  be  adjusted. 

The  structural  matrices  for  the  rotating  system  were  obtained  by  examining  the 
expression  for  the  total  system  energy.  This  provides  a  differential  stiffness  matrix  that 
models  centrifugal  force  and  a  provides  velocity -dependent  Coriolis  matrix.  By  using  a 
high-level  programming  language  (Direct  Matrix  Abstraction  Programming)  a  modal 
analysis  solution  sequence  was  modified  to  account  for  rotational  effects  in  free  vibration. 
Finite  element  models  were  then  created  for  a  typical  compressor  blade  in  a  modern  jet 
engine  and  for  a  cantilever  beam  rotating  about  the  vertical  axis. 1  Using  these  models,  the 
effects  of  rotation  as  simulated  by  the  finite  element  method  were  verified  against 
theoretical  results. 

The  optimal  redesign  was  done  by  deriving  complex  nonlinear  inverse  perturbation 
equations  for  the  problem  involving  both  magnitude  and  phase  components.  The 
perturbation  problem  is  solved  by  using  nonlinear  mathematical  programming.  In  order  to 
optimall^redesign^^n^  uses  an  underdetermined  system,  i.e.,  the  feasible  design  must  not 
be  unique.  This  allows  the  application  of  an  objective  function,  such  as  minimum 
structural  weight  or  minimum  change  from  the  baseline  design.  Constraints,  such  as  those 
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CHAPTER  1 


INTRODUCTION 

The  desire  to  minimize  weight  while  meeting  design  requirements  is  a  key  concept 
in  optimal  structural  design.  When  a  particular  stuctural  component  is  being  examined, 
the  engineer  must  consider  many  aspects  of  structural  behavior,  such  as  dynamic 
response.  These  parameters  must  be  bounded  in  the  design  phase  while  still  meeting  the 
overall  goal  of  minimizing  weight. 

There  are  numerous  occasions  where  the  design  of  an  existing  structure  is  found 
wanting  and  this  creates  the  need  for  structural  redesign.  For  example,  the  mission  of  a 
particular  aircraft  that  has  been  used  for  many  years  may  change.  This  happened  with 
the  T-38  trainer  when  it  was  adapted  for  use  in  lead-in  fighter  training.  The  aircraft 
structure  may  then  be  subjected  to  loads  that  could  not  be  forseen  in  the  original  design. 
Fatigue  cracks  may  develop  due  to  changes  in  the  load  spectrum.  Dynamic  effects  may 
result  due  to  excitation  of  modal  frequencies  that  would  have  remained  stable  under  the 
prior  usage. 

In  any  structural  redesign  problem,  there  are  several  possible  candidate  designs  that 
would  meet  the  new  criteria.  Structural  optimization  can  be  used  in  the  redesign  problem 
to  find  the  best  possible  configuration.  The  objective  function  for  redesign,  however,  may 
be  quite  different  from  the  objective  function  in  the  original  design  of  a  part.  It  may  not  be 
suitable  to  concern  oneself  only  with  weight.  If  the  dimensions  of  an  existing  part  are 
altered  too  much,  it  may  r'  longer  fit  in  the  aircraft.  Existing  fastener  locations  may  no 
longer  be  suitable.  The  balance  of  the  entire  aircraft  may  be  thrown  off.  Therefore,  it 
may  be  necessary  to  obtain  the  optimal  design  that  requires  the  least  possible  change  from 


the  original  design  in  terms  of  some  dimension,  such  as  thickness.  This  is  a  minimum 
design  change  criterion  for  optimization. 

For  given  values  of  design  parameters,  analytical  methods  may  be  used  to 
uniquely  determine  the  response.  The  response  values  can  be  compared  against  the 
desired  response.  A  search  in  the  design  space  can  be  undertaken  for  a  new  design,  and 
the  response  and  constraints  can  be  reanalyzed.  This  approach  is  the  simplest 
conceptually,  and  if  the  number  of  design  variables  and  constraints  is  small,  is  the  most 
advantageous.  However,  most  problems  in  structural  design  involve  a  large  number  of 
design  variables  and/or  a  large  number  of  constraints.  This  makes  a  random  search 
approach  unfeasible.  In  addition,  a  closed  form  analytical  solution  usually  cannot  be 
obtained  for  a  large  problem.  Large,  general  purpose  finite  element  programs  are  used  to 
determine  the  structural  response.  Since  these  programs  are  expensive  to  run  on  a 
computer,  it  is  desirable  to  avoid  large  numbers  of  trial  solutions.  Therefore,  some  form  of 
optimization  procedure  coupled  with  finite  element  analysis  must  be  used. 

Some  structures,  because  of  their  application,  are  subject  to  body  forces  in  addition 
to  the  loadings  caused  by  boundary  forces.  Rotating  bodies  experience  centrifugal  effects. 
These  effects  are  not  forces  at  all  if  viewed  from  a  nonrotating  frame  of  reference. 
However,  in  the  rotating  frame,  the  centrifugal  effects  may  be  viewed  as  a  reverse- 
effective  pseudo-force.  If  the  rotational  speed  is  small  enough,  these  forces  can  be 
neglected.  In  high  speed  applications,  however,  the  body  forces  must  be  taken  into 
account.  Centrifugal  force  in  particular  can  result  in  stiffening.  This  stiffening  is 
particularly  obvious  for  rotating  helicopter  blades.  When  the  blades  are  at  rest,  they  sag 
under  their  own  weight.  As  the  blades  speed  up,  they  stiffen,  to  the  point  where  they  can 
bear  not  only  their  own  weight,  but  the  weight  of  the  entire  aircraft.  A  similar  effect 
occurs  in  aircraft  gas  turbine  engines.  The  rotating  blades  in  the  high  speed  compressor 
and  turbine  are  subject  to  body  forces.  With  the  advent  of  higher  speed  blading,  it  is 
necessary  to  include  the  effect  of  these  body  forces  in  design.  This  introduces  a 
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nonlinearity  in  the  finite  element  analysis  and  optimization  scheme. 

Coriolis  force  is  another  body  “force”  that  a  rotating  body  may  experience  when 
viewed  in  a  rotating  coordinate  system  but  is  not  a  force  at  all  if  viewed  in  an  absolute 
frame  of  reference.  This  force,  also  known  as  gyroscopic  force,  couples  motion  in  one  plane 
with  motion  in  another  plane.  Gyroscopic  effects  are  velocity  dependent;  that  is,  the 
greater  the  velocity  in  one  plane,  the  greater  the  effect  will  be  in  the  other  plane.  In 
systems  which  permit  large  amounts  of  out-of-plane  motion,  such  as  in  pretwisted  blades, 
the  Coriolis  forces  will  become  great  at  high  speeds,  and  thus  will  affect  optimal  redesgin. 
In  many  other  systems,  Coriolis  forces  may  be  ignored.  The  effect  of  the  Coriolis  terms  on 
structural  dynamics  is  documented  by  Greenwood  (1977)  and  Meirovitch  (1980). 

The  effect  of  rotation  on  the  dynamic  characteristics  of  a  body  has  been  well- 
documented  by  Washizu  (1982  ).  Simplified  methods  for  the  solution  of  the  problem  of  the 
rotating  beam  have  been  presented  by  White  and  Malatino  (1975),  McDaniel  and  Murthy 
(1977),  and  Giurgiutiu  and  Stafford  (1977).  Isakson  and  Eisley  (1964)  related  the 
centrifugal  effect  to  a  stiffening  parameter,  the  Southwell  coefficient.  Hodges  and 
Rutkowski  (1981)  derived  a  finite  element  solution  to  the  problem  of  a  rotating  beam  by 
creating  variable  order  shape  functions.  Queau  and  Trompette  (1981)  used  parameter 
optimization  in  a  beam  analysis  of  a  blade.  Cross-sectional  properties  were  found  under  a 
least  weight  objective  and  with  frequency  and  stress  constraints.  Kim,  Anderson,  and 
Sandtrom  (1983)  presented  a  nonlinear  inverse  perturbation  scheme  for  automated 
redesign  of  modal  characteristics.  Hoff,  et.  all  (1984)  used  the  perturbation  technique  to 
develop  a  predictor-corrector  approach,  which  is  useful  in  the  present  study. 

The  advance  of  numerical  techniques  of  structural  optimization  has  been  one  of  the 
most  important  developments  in  the  field  of  engineering  design.  These  methods  have 
permitted  the  engineer  to  select  the  best  possible  configuration  from  several  candidate 
designs.  Structural  optimization  can  also  be  used  to  modify  an  existing  design  in  order  to 
adapt  existing  aircraft  to  new  requirements.  One  critical  area  where  optimal  structural 
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redesign  will  be  of  great  importance  is  in  the  reanalysis  of  aircraft  gas  turbine  blading.  As 
the  demand  for  increased  performance  continues,  blades  may  have  to  be  redesigned  to 
counteract  the  effect  of  loading  changes.  New  procedures  will  have  to  be  created  to  do  this 
efficiently. 

The  objective  of  this  study  is  to  first  create  a  finite  element  solution  for  the  stiffening 
effect  of  centrifugal  force  in  order  to  calculate  the  mode  shapes  and  frequencies  of  a 
structure  composed  of  several  different  types  of  elements.  Then,  nonlinear  inverse 
perturbation  will  be  used  to  account  for  the  changing  nature  of  the  stiffening  effect  as  the 
design  process  progresses.  This  scheme  measures  design  changes  as  perturbations  of  the 
original  structure.  The  method  will  be  applied  to  a  problem  involving  a  typical  curved 
compressor  blade.  Finally,  the  effects  of  Coriolis  forces  on  the  optimal  redesign  process 
will  be  studied. 

This  dissertation  begins  with  a  review  of  rotating  systems,  structural  optimization, 
and  the  optimization  of  rotating  systems.  The  effect  of  rotation  on  frequencies  and  mode 
shapes  is  discussed  in  Chapter  3.  A  finite  element  procedure  to  account  for  static  loads  in 
dynamic  response  is  also  described  in  detail.  The  effect  of  Coriolis  forces  is  also 
considered. 

In  Chapter  4,  the  rotational  characteristics  of  rotating  systems  are  examined.  The 
stiffening  effect  of  centrifugal  forces  is  examined  both  by  a  finite  element  approach  and  by 
a  classical  theoretical  formulation. 

In  Chapter  5,  a  nonlinear  inverse  perturbation  scheme  including  centrifugal  and 
Coriolis  effects  is  developed.  The  application  of  Sequential  Unconstrained  Minimization 
Technique  (SUMT)  to  the  problem  in  optimization  is  described. 

The  predictor-corrector  technique  for  the  optimization  of  rotating  systems  is  derived 
in  Chapter  6.  The  method  is  applied  to  several  simple  examples. 

Chapter  7  applies  the  procedure  to  a  large  problem.  A  fiat  compressor  blade  is  used 
as  an  example  to  check  the  method.  The  effect  of  Coriolis  forces  on  the  optimization 
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problem  is  examined  in  a  rotating  beam. 
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CHAPTER  2 

LITERATURE  REVIEW 

Literature  Survey  on  Rotation 

The  characteristics  of  vibrating  systems  were  treated  mathematically  first  by 
Rayleigh  (1898).  The  basic  equations  of  motion  of  a  rotating  beam  of  constant  section  with 
axis  of  rotation  perpendicular  to  the  length  of  the  beam  were  updated  by  Boyce.  Di  Prima, 


and  Handelman  (1954).  Nagaraj  and  Shanthakumar  (1975)  solved  the  problem  with 
Galerkin’s  method.  Tomar  and  Dhole  (1976)  expanded  the  solution  to  encompass 
pretwisted  beams,  solving  the  equations  by  a  conventional  Rayleigh-Ritz  approach.  Putter 
and  Manor  (1978)  devised  a  finite  element  model  of  a  rotating  tapered  beam  that 
incorporated  centrifugal  effects.  Hoa  (1979)  added  a  tip  mass  to  the  finite  element 
procedure. 

The  basic  equations  of  motion  for  a  rotating  cantilever  blade  have  been  extensively 
investigated.  Lo  and  Renbarger  (1951)  analyzed  rotating  beams.  Houbolt  and  Brooks 
(1958)  showed  that  centrifugal  force  increases  the  first  bending  frequency  for  a  rotating 
blade.  Carnegie  (1959)  and  Carnegie  (1967),  using  energy  methods,  expanded  this  work. 
Isakson  and  Eisley  (1964)  determined  the  natural  frequencies  of  twisted  blades  for  both 
rotating  and  nonrotating  systems.  This  subject  was  also  studied  by  Carnegie  and  Dawson 
for  straight  blades  (1969)  and  for  pretwisted  blades  (1971).  Various  solution  methods  for 
these  equations  have  been  proposed.  Wadsworth  and  Wilde  (1967)  solved  the  problem 
using  a  pair  of  simultaneous  differential  equations  and  a  Runge-Kutta  approach.  Rao  and 
Carnegie  (1970)  used  a  Ritz  approach.  A  transfer  matrix  approach  was  outlined  by 
McDaniel  and  Murthy  (1977).  A  mixed  variational  approach  was  outlined  by  Lang  and 
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Nemat-Nasser  (1979).  An  integrating  matrix  finite  difference  procedure  was  described  by 
Hunter  (1970).  An  improved  method  for  the  analysis  of  pretwisted  airfoil  blades  was 
presented  by  Sybrahmanyam  and  Kaza  (1984). 

Coriolis  effects  have  also  been  discussed  by  several  authors.  Hunter  (1970) 
dismisses  these  effects  as  being  insignificant  in  comparison  with  centrifugal  forces; 
however.  Subrahmanyam  and  Kaza  (1985)  determined  that  Coriolis  forces  are  an 
important  influence  at  high  rotational  speeds  and  for  thick  beams.  They  also  stated  that 
Coriolis  forces  may  be  important  even  at  low  speeds  under  certain  geometries.  This 
analysis  was  expanded  by  Subrahmanyam  et  al  (1987)  to  include  pretwisted  blades. 
Anarsi  (1986)  also  includes  the  effects  of  shear  deformation  and  rotary'  inertia.  Sisto  et  al 
(1983)  analyzes  the  stability  of  a  rotating  blade  by  the  use  of  Floquet  theory. 

Finite  element  solutions  to  the  problem  of  free  vibration  of  rotating  beams  was 
examined  by  Hodges  and  Rutkowski  (1981).  Thomas  and  Subuncu  (1979),  using  a  finite 
element  method,  solved  the  vibrational  characteristics  of  rotating  pretwisted  asymmetric 
cross-sectional  blading.  Subuncu  (1985)  produced  a  finite  element  solution  for  blades  with 
a  nonlinear  angular  pretwist;  however  his  method  did  not  involve  rotational  effects. 
Dokainish  and  Rawtani  (1971)  obtained  a  finite  element  solution  for  a  rotating  cantilever 
plate.  They  included  centrifugal  forces  but  neglected  Coriolis  effects.  A  textbook  on  the 
topic  was  published  by  MacNeal  (1973). 

Literature  Survey  in  Optimization 

The  problem  of  optimal  structural  design  has  been  studied  for  many'  years.  At  first, 
a  trial-and-error  approach  was  used,  with  initial  choice  based  on  experience  (Sheu  and 
Praeger,  1968).  However,  methods  of  optimization  now  exist  that  obtain  the  desired 
structural  configuration  more  efficiently.  These  approaches  may  be  categorized  as  either 
direct  or  indirect  methods  (Kiusalaas,  1972).  In  indirect  methods,  an  optimality  criteria  is 
used.  Direct  methods  are  based  on  mathematical  programming. 

The  use  of  optimality  criteria  lies  at  the  heart  of  the  indirect  approach  to  structural 


optimization.  Important  developments  in  this  procedure  have  been  made  by  Prager  and 
Taylor  (1967)  and  Taylor  (1968).  These  methods  involve  the  establishment  of  a  criterion 
that  defines  the  optimum,  such  as  uniform  strain  energy  density,  and  the  creation  of  an 
iterative  procedure  for  obtaining  the  design  (Venkayva,  Khot,  and  Berke,  1973).  The 
Kuhn-Tucker  conditions  are  the  first  order  necessary  conditions  for  obtaining  an  acceptable 
design  (Luenberger,  1984). 

The  direct,  or  mathematical  programming,  methods  use  numerical  search  techniques 
to  solve  the  problem  of  optimal  structural  design.  Moses  (1964)  presented  a  method  for 
finding  the  optimum  design  by  using  sequential  linear  programming.  A  more  general 
approach  is  to  minimize  the  objective  function  as  an  unconstrained  function,  but  to  impose 
a  penalty  in  order  to  limit  constraint  violations.  This  requires  that  several  unconstrained 
minimization  problems  be  solved.  One  such  method,  called  the  Sequential  Unconstrained 
Minimization  Technique  (SUMT)  was  extensively  described  by  Fiacco  and  McCormick 
(1968). 

Derivatives  of  structural  response  with  respect  to  the  design  variables  can  be 
obtained  and  used  in  gradient  projection  methods  (Fletcher  and  Reeves,  1964).  These 
derivatives  can  be  calculated  by  a  design  space  method  (Fox,  1965).  Kavanaugh  (1972) 
used  dynamic  relaxation  to  devise  an  approximate  algorithm  for  uncoupling  the 
optimization  problem. 

The  application  of  structural  redesign  to  dynamic  problems  has  been  extensively 
investigated.  Rayleigh  (1898)  provided  the  first  solution  for  designing  a  dynamic  system  to 
meet  specific  frequency  requirements.  Turner  (1967)  used  the  equations  of  motion  for  a 
large  deflection  vibration  problem  as  equality  constraints  and  employed  a  Lagrange 
multiplier  method.  He  used  a  finite  element  idealization  and  an  iterative  scheme  to  solve 
the  nonlinear  equations.  Taylor  (1967)  proposed  an  alternate  method.  He  developed  a 
functional  that  related  the  system  energy  to  the  eigenvalue  problem,  and  used  this 
functional  to  obtain  the  optimum  design.  Taylor  (1968)  added  a  condition  that  the  cross- 
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sectional  area  not  be  reduced  below  a  specified  minimum  during  the  redesign  process.  This 
was  accomplished  by  an  inequality  constraint  and  Lagrange  multipliers.  Sheu  (1968) 
applied  the  approach  of  Turner  and  Taylor  to  a  one-dimensional  structure  and  added 
segment  boundaries  and  specific  stiffnesses  as  design  requirements.  Taylor  (1969)  applied 
an  energy  formulation  to  the  problem  of  truss  design.  Sippel  and  Warner  (1973)  devised  a 
mimimum  mass  optimality  criterion  and  compared  the  results  of  continuous  optimum 
design  to  piecewise  uniform  optimum  design. 

McCart,  Haug,  and  Streeter  (1970)  used  a  steepest  descent  approach  to  obtain  the 
optimal  design  for  a  three-bar  frame.  Rubin  (1970)  satisfied  a  frequency  constraint  and 
found  the  minimum  structural  weight  through  a  step-wise  procedure  that  defined  two 
“modes  of  travel”.  The  first  mode  changed  the  frequency  to  obtain  a  desired  value.  The 
second  one  minimized  structural  weight  while  keeping  the  frequency  constant.  The 
problem  with  this  approach  is  that  while  it  may  obtain  a  design  it  may  not  necessarily  find 
the  best  design.  Armand  (1971)  used  methods  from  optimal  control  theory  to  present  a 
different  approach  to  the  optimal  design  of  two-dimensional  structures. 

Taylor  (1977)  obtained  a  procedure  that  upgrades  a  computer  model  of  a  structure 
so  that  the  calculated  normal  mode  frequencies  more  nearly  match  those  determined  from 
experiment.  This  method  was  based  on  a  Taylor  series  type  expansion.  Only  first  order 
terms  were  retained. 

If  a  baseline  structure  exists  but  it  is  found  that  the  response  of  the  structure  is 
unacceptible  and  modification  is  necessary,  perturbation  techniques  may  be  employed  to 
obtain  the  desired  values.  Stetson  (1975)  introduced  small  changes  in  mass  and  in  the 
stiffness  moduli  of  a  structure.  He  then  obtained  a  first  order  perturbation  method  that 
obtained  the  mode  shapes  for  the  perturbed  structure.  He  introduced  the  concept  of 
“admixture  coefficients”  that  expressed  the  mode  shapes  of  the  perturbed  structure  in 
terms  of  combinations  of  the  baseline  mode  shapes.  Stetson  and  Palma  (1976),  Stetson  et 
al  (1978),  and  Stetson  and  Harrison  (1981)  expanded  this  technique  to  encompass  a  finite 
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element  structural  formulation  and  applied  it  to  several  problems.  Sands trom  and 
Anderson  (1982)  related  Stetson’s  admixture  coefficients  to  physical  changes  in  the  finite 
element  model.  Kim  et  al  (1983)  obtained  a  complete  nonlinear  inverse  perturbation 
technique  using  the  equations  of  dynamic  equilibrium.  A  SUMT  penalty  function  method 
was  used  where  the  objective  function  was  minimum  weight  and  the  penalty  term  involved 
a  normalized  set  of  residual  force  vectors. 

One  major  problem  of  the  nonlinear  inverse  perturbation  method  is  that  for  a  large 
problem,  the  number  of  calculations  required  become  excessive.  For  that  reason,  Kim  and 
Anderson  (1984)  and  Kim  (1985)  used  generalized  dynamic  reduction  to  transform  the 
problem  into  a  small  sized  subspace.  Hoff  et  al  (1984)  overcame  the  difficulties  in  applying 
the  nonlinear  inverse  perturbation  method  by  using  an  incremental  predictor-corrector 
technique.  In  the  predictor  phase,  element  changes  necessary'  to  enforce  the  desired  mode 
shape  and  frequency  changes  are  obtained  through  a  first  order  solution  of  the  dynamic 
equations.  In  the  corrector,  approximate  eigenvectors  are  obtained  for  the  objective 
system,  which  are  then  used  to  correct  the  elemental  changes.  This  method  was  later 
expanded  by  Hoff  (1985). 

Literature  Survey  in  the  Optimization  of  Rotating  Systems 

Olhoff  and  Parbery  (1982)  analyzed  the  design  of  vibrating  beams  and  rotating 
shafts  with  lumped  mass  at  the  tip  and  at  midspan.  Their  method  incorporated 
centrifugal,  but  not  Coriolis,  effects.  Kengtung  and  Gu  (1984)  also  studied  the  problem  of 
rotating  blades  without  Coriolis  effects  for  small  (ten  element)  problems  using  beam 
elements.  Kounadis  (1985)  analyzed  the  effect  of  axial  inertia  on  the  bending  frequencies 
of  a  frame  structure. 

Bennett  (1983)  examined  the  the  application  of  optimization  methods  to  problems  in 
the  design  of  helicopter  rotors.  His  optimization  scheme  involved  a  linear  design 
sensitivity  approach.  His  method  did  not  account  for  changes  in  blade  weight  and  inertia 
properties  since  these  affect  centrifugal  forces  in  a  nonlinear  manner.  Therefore,  he 
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required  that  blade  weight  and  inertia  remain  fixed  during  the  process  of  optimal  design. 

Queau  and  Trompette  (1981)  applied  changes  in  inertia  properties  during  the 
redesign  process  to  determine  changes  in  centrifugal  stiffening  affecting  optimization. 
Their  method  also  involved  linear  design  sensitivites.  It  also  did  not  update  the 
eigenvectors  during  the  optimization  process,  requiring  a  large  number  of  calculations. 
Only  beam  elements  were  used.  Their  method  of  optimization  can  lead  to  a  local 
minimum.  Their  procedure  cannot  be  used  on  a  general  class  of  problems,  particularly 
plate-like  bodies. 


CHAPTER  3 


BACKROUND  OF  ROTATIONAL  EFFECTS 

In  this  chapter,  the  equations  of  motion  are  discussed  for  small  motions  of  a  rotating 
discrete  system.  The  terms  in  the  equations  that  represent  rotational  effects  are  shown  in 
a  structural  mechanics  representation.  The  corresponding  eigenvalue  problem  and  its 
methods  of  solution  are  shown.  Finally,  a  physical  example  will  be  used  to  derive  the 
Coriolis  matrix  for  a  rotating  bar. 

Equations  of  Motion  of  a  Rotating  Discrete  Structure 
A  system  consisting  of  n  degrees  of  freedom  has  its  motion  fully  described  by  n 
generalized  coordinates  qt(t),  where  i=  l,2,...,n.  The  kinetic  energy  of  the  system,  T,  may 
be  defined  as  (Meirovitch,  1980): 

T  =  T0  +Tj  +  T2  (3.1) 

where 

t2=|e  i>.M  <s-2) 

1=1  j= i 

is  a  quadratic  function  of  the  generalized  velocities  q,(t). 
and 

T,  =  £f,q,  (3-31 

i=  l 

is  a  linear  function  of  the  generalized  velocities.  The  term  T0  contains  no  generalized 
velocities:  however,  the  function  T0  and  the  coefficients  and  f,  will  depend  on  the 
generalized  coordinates  q,.  The  term  T0  gives  rise  to  terms  involving  centrifugal  force  and 
behaves  as  a  potential  energy.  The  Tj  term  produces  Coriolis  forces. 
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In  addition  to  centrifugal  and  Coriolis  forces,  there  may  be  other  forces  that  act  on 
the  system.  These  include  forces  from  a  potential  energy  function  written  as  a  function  of 
the  generalized  coordinates  alone,  with  no  time  derivatives: 

V  =  V(qs)  (3.4: 

Such  forces  can  be  elastic  forces  or  gravitational  forces.  These  elastic  forces  generate  the 
conventional  small-displacement  system  stiffness. 

The  final  set  of  forces  are  those  that  fall  into  no  particular  set,  and  shall  be 
represented  by  Qi(  the  generalized  forces.  These  forces  may  depend  on  time  but  not  on 
displacement  or  velocities.  Friction  forces  are  an  example  of  generalized  forces. 

The  equations  of  motion  of  the  system  may  be  found  through  the  use  of  the  classic 
differential  form  of  the  Lagrange  equation: 

dL 

- =Qi  (3.5: 

dt  \dqj  dqi 

where  the  Langrangian  L  is  defined  by: 

L  =  T  -  V  (3.6: 

The  system  of  equations  represented  by  Equation  (3.5)  comprises  a  set  of  n 
nonhomogeneous  nonlinear  ordinary  differential  equations  of  second  order.  General 
solutions  of  this  type  of  equations  do  not  exist.  Under  certain  conditions  it  is  possible  to 
make  simplifying  assumptions  that  will  permit  the  linearization  of  the  equations.  One  such 
assumption  is  that  motion  will  be  restricted  to  small  motions  in  the  i.  igborhood  of  the 
static  equilibrium  condition.  qi0  =  0,  for  i  =  1.2.. ..,n.  A  coordinate  transformation  can  be 
done  to  translate  the  origin  so  as  to  make  it  coincide  with  the  static  equilibrium  point. 

Thus,  only  motions  in  the  neighborhood  of  the  trivial  solution  ql0  =  Qio  =  ®  (>=  L2 . n)  will 

be  considered. 

Using  these  assumptions,  it  can  be  shown  (Meirovitch,  1980i  that: 
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mij  =  mji=- 
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(3.7) 


dqj&ii 


where  the  terms  my  are  called  the  mass  or  inertia  coefficients. 


It  can  also  be  shown  that: 

where  the  terms  fy  are  constant  coefficients.  This  implies  that: 


f«a  — 


(3.8) 


T»*EIW 

I=lj=l 


(3.9) 


Let  U  be  defined  by: 


U  =  V  -  Tn 


(3.10) 


where  U  represents  a  dynamic  potential.  Furthermore,  the  stiffness  coefficients  ky  can  be 
shown  to  be: 


32U 


^ij  ^ji 


2 

dU 


(3. IF 


dqjdqi 


As  discussed  below,  the  centrifugal  effects  can  be  treated  as  a  static  preload  on  the 
system.  This  implies  that  centrifugal  forces  will  be  included  in  the  strain  energy7 
component  V  of  the  potential  energy  U.  rather  than  as  T0  kinetic  energy.  A  full  derivation 
of  the  tangent  stiffness  matrix  including  centrifugal  effects  is  shown  in  Chapter  4. 

The  explicit  equations  of  motion  can  be  derived  and  expressed  as: 


^[myqj  +  byqj  +  kyqj]  =  Q,  (i=l,2,...n) 
j=i 


(3.121 


Note  that  the  mass  and  stiffness  coefficients  are  symmetric  while  the  Coriolis  coefficients, 
by,  are  skew  symmetric  such  that: 
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bjj  =  fjj  -  fy,  (ij=l,2,...,n) 


(3.13) 


In  matrix  form,  Equation  (3.12)  may  be  written  as: 


PVI]{q}  +  [B]{q>  +  [K]{q}  =  {Q} 


(3.14) 


The  symmetry  of  the  coefficients  and  ky  implies  that: 


—  nwiT 


[M]  =  [M] 


(3.15) 


—  nriT 


[K]  =  [K] 


(3.16) 


while  the  skew  symmetry  of  the  coefficients  by  implies  that: 
[B]  =  -[B]t 


(3.17) 


Another  way  to  deal  with  the  differences  between  centrifugal  forces  and  Coriolis 


effects  is  as  follows.  In  structural  analysis,  centrifugal  force,  which  is  dependent  only  on 


mass,  position,  and  rotational  speed,  can  be  considered  as  a  static  preload  on  the  system. 


As  such,  its  effects  can  be  included  in  the  initial  stress,  or  differential  stiffness,  matrix. 


Coriolis  effects,  however,  are  dependent  on  velocity,  which  makes  them  dependent  on  the 


first  derivative  of  the  eigenvector.  Coriolis  terms,  unlike  structural  damping,  are  energy- 


conserving.  Just  as  the  structural  matrix  derived  from  centrifugal  forces  can  be  absorbed 


into  the  stiffness  matrix,  the  matrix  derived  from  Coriolis  effects  could  be  absorbed  into 


the  system  damping  matrix  if  there  is  one. 


Eigenvalue  Problem  for  Conservative  Coriolis  Systems 


For  the  free  vibration  problem,  the  damping  and  the  generalized  forces  Q,  are  taken 


to  be  zero.  Therefore,  Equation  (3.12)  may  be  written  by: 


[M]{q]  +  [B]{q}  +  [K]{q}  =  {0} 


(3.18) 


MacNeal  (1973)  showed  that  the  system  described  in  Equation  (3.18)  is  energy- 


conserving  if  the  mass  and  stiffness  matrices  are  real,  symmetric,  positive  definite  and  if 


the  Coriolis  matrix  is  real  and  skew  symmetric.  The  solution  to  Equation  (3.18)  will  be  of 


the  form: 


{q(t)}  =  epl{q} 


(3.19) 


where  p  is  a  constant  complex  scalar  and  {q}  is  a  constant  complex  vector.  Equation  (3.19) 
is  introduced  into  Equation  (3.17)  and  ept  is  canceled  to  obtain  the  eigenvalue  problem: 

(p2[M]  +  p[B]  +  (K]){q}  =  0  (3.20) 

The  characteristic  equation  for  this  system  is  obtained  by  setting  the  determinant  of 
the  coefficients  equal  to  zero: 

|  p2[M]  +  p[B] + [K]  |  =  0  (3.21) 

This  equation  gives  a  polynomial  of  degree  2n  in  p.  Greenwood  (1977)  states  that  due  to 
the  symmetry  of  the  mass  and  stiffness  matrices  and  the  skew  symmetry  of  the  Coriolis 
matrix,  all  of  the  odd  powers  of  p  are  absent  from  the  characteristic  equation.  Therefore, 
a  polynomial  of  nth  degree  in  p2  is  generated.  This  would  also  be  the  case  if  the  Coriolis 
terms  were  absent.  Therefore,  the  eigenvalues  will  consist  of  n  pure  imaginary  conjugate 
pairs,  pr=  ±iwr,  where  r=  l,2,...,n.  The  eigenvectors  will  also  occur  in  complex  conjugate 
pairs,  {q}r  =  {y}r  +  i{z},  {q}r  =  {y}r  -  i{z}r,  where  {y}r  is  the  real  part  and  i{z}r  is  the 
imaginary  part  of  the  eigenvector  {q]r.  This  implies  that  as  a  result  of  the  Coriolis  effects, 
the  amplitude  ratios  will  not,  in  general,  be  real.  Therefore,  the  components  of  an 
eigenvector  pair  for  a  given  eigenvalue  pair  will  oscillate  at  the  same  frequency  but  not  in 
phase. 

Meirovitch  (1974)  provides  an  alternate  solution  to  the  system  of  Equation  (3.  IS) 
tha.  reduces  it  to  a  standard  form  that  is  similar  to  the  system  without  Coriolis  forces.  If 
the  matrices  in  Equation  (3.18)  are  of  order  n,  then  one  can  introduce  a  2n-dimensional 
state  vector: 

{x}T  =  [{q}|{q)]  (3-22) 

Therefore,  Equation  (3.18)  can  be  rewritten  as: 

[M*]{x}  +  [B*]{x}  =  {0}  (3.23) 


where 


are  real  nonsingular  matrices  of  order  2n. 

Though  this  method  reduces  the  S3'Stem  of  equations  involving  Coriolis  components 
to  a  more  conventional  form,  the  arrays  created  are  twice  as  large.  This  will  force  a 
heavy  penalty  in  storage  requirements  in  numerical  implementation.  In  actual  application, 
the  direct  solution  of  the  complex  problem  posed  in  Equation  (3.18)  is  more  efficient  than 
the  solution  of  the  modified  system  of  Equation  (3.23).  Therefore,  this  method  posed  by 
Meirovitch  will  not  be  used. 


Derivation  of  Coriolis  Matrix  for  a  Bar  Element 
The  Coriolis  matrix  for  a  rotating  bar  element  will  now  be  derived,  using  a  lumped 
mass  approach.  Figure  1  shows  a  typical  rotating  bar.  The  entire  mass  of  the  bar  will  be 
divided  into  two  equal  lumped  masses  at  each  of  the  two  nodes.  Each  mass  m  will  have 
one-half  of  the  mass  of  the  total  mass  of  the  bar.  The  bar  rotates  around  the  global  z-axis 
at  a  rotational  speed  of  f!  hz. 

The  position  vector  to  the  mass  at  the  first  node,  {rj}.  is  given  by: 


{r1}  =  x1i  +  yij  +  z,k 


(3.24) 


where  Xj,  ylt  and  Zj  are  the  coordinates  of  node  1  in  the  global  system  and  i.  j,  and  k  are 
the  unit  direction  vectors  for  the  x,  y,  and  z  axis  respectively.  Similarly,  the  position 
vector  to  the  mass  at  the  second  node  {r2}  is  given  by: 


{r2}  =  x2r+yj+z2k 


(3.25) 


The  absolute  time  derivative  of  a  vector,  {A},  in  a  rotating  system  is  given  by: 


Figure  1.  Rotating  Bar  -  Lumped  Mass 


{A}0  =  {A}r  +  {n}x{A}  (3.26) 

where  {A}0  is  the  absolute  time  derivative  of  {A},  {A}r  is  the  time  derivative  of  {A}  with 
respect  to  the  rotating  system,  and  {fi}  is  the  rotation  vector. 

It  is  assumed  that  the  rotation  will  be  purely  about  the  z  axis.  Thus  the  rotation 
vector  is  given  by  Ok.  Therefore,  the  absolute  time  derivative  of  the  position  vector  to 


node  i,  also  called  the  velocity  {v}^ 

{v)i  =  (Xj  -  fty,)i + (y,  +  +Zjk 


(3.27) 


To  find  the  acceleration  of  the  mass  at  node  i,  Equation  (3.26)  is  applied  to  Equation 


T 


(3.27)  to  obtain: 


{a}i  =  (xi-2fiy1-fi2xi)i  +(yj  +  2f2Xj-  +z;k 


(3.28) 


The  second  time  derivitive  terms  in  the  above  equation  represent  the  acceleration  of 


the  mass  i  relative  to  the  rotating  frame.  The  terms  dependent  on  Q2  are  the  centrifugal 


acceleration  components.  Notice  that  they  are  dependent  only  on  the  position  of  the  mass, 


with  x-displacement  giving  the  dependency  of  centrifugal  acceleration  in  the  x-direction  and 


y-displacement  giving  the  dependency  of  centrifugal  acceleration  in  the  y-direction.  The 


components  of  Equation  (3.26)  dependent  on  fi  are  the  Coriolis  terms.  It  is  important  to 


note  that  Coriolis  acceleration  in  the  x-direction  depends  on  the  velocity  of  the  mass  i  in  the 


y-direction  and  that  Coriolis  acceleation  in  the  y-direction  depends  on  velocity  in  the  x- 


direction. 


Relative,  centrifugal,  and  Coriolis  acceleration  can  each  be  considered  to  be  reversed- 


effective  forces  or  D’Alembert  forces.  As  shown  above,  the  reversed-effective  forces  due  to 


relative  acceleration  and  centrifugal  acceleration  in  the  x-direction  will  be  dependent  upon 


acceleration  or  displacement  in  the  x-direction,  respectively,  with  similar  forces  in  the  y- 


direction.  However,  the  reversed-effective  Coriolis  force  that  points  in  the  x-direction 


depends  on  y-direction  velocity,  and  vice-versa. 


The  kinetic  energy  of  the  mass  at  the  ith  node  is  given  by: 


T1=|miM1'{v}1 


(3.29) 


T‘=|m[(x1-nyl)2  + (y^^)2  +z2] 

Equation  (3.30)  can  be  broken  down  into  its  components,  T0  through  T2.  to  provide 


(3.30) 


the  desired  structural  matrices.  The  kinetic  energy  components  are: 


T>2  =  i{q}T[Ml{q} 


(3.31) 


T\=i{q}T[F]{q} 


(3.32) 
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where 


[M]  = 

i 


0 

0 


0 

mi 

0 


0 

0 

mi 


I 

!  [F]= 


0 

n\jQ 

0 


—  mjfi 
0 
0 


0 

0 

0 


i 

'  [S]  = 


mjO2  0 
0  nijO2 

0  0 


0 

0 
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To  form  the  Coriolis  matrix  from  [F],  one  remembers  from  Equation  (3.13)  that: 


[B]  =  [F]t  -  [F]  (3.34) 

where  [F]  is  the  matrix  of  the  f^  factors. 

In  order  for  one  to  form  element  matrices  from  the  matrices  given  in  Equations 
(3.31)  through  (3.34),  the  components  from  the  other  node  must  be  added  in.  Thus  [M] 
will  become  the  element  mass  matrix,  [B]  will  become  the  element  Coriolis  matrix,  and  [S] 
will  become  the  element  stiffness  matrix  from  the  applied  load  due  to  the  rotation.  This 
last  matrix,  however,  differs  from  the  differential  stiffness  matrix.  Differential  stiffness 
results  from  examining  small  nonlinear  terms  in  the  derivation  of  stiffness  to  form  the 
tangential  stiffness  matrix  for  the  element.  These  nonlinearities  result  from  applied 
loading,  such  as  centrifugal  force. 

The  matrix  derived  from  the  T0  component  of  kinetic  energy  has  been  previously 
used  to  model  the  effect  of  centrifugal  loading  on  the  eigenfrequencies  of  rotating 
structures  (Trompette  and  Lalanne,  1974,  Olhoff  and  Parbery,  1982,  Subrahmanyam  and 
Kaza,  1984,  Queu  and  Trompette,  1981).  In  the  current  analysis,  however,  this  approach 
will  not  be  used.  Instead,  the  forces  generated  by  the  centrifugal  effect  will  be  considered 
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as  an  applied  load  (Hoa,  1978).  Optimization  can  then  be  done  in  the  presence  of  the 


applied  loads  (Garstecki,  1984).  A  potential  energy  expression  will  then  be  derived  from 


these  forces  (MacNeal,  1972).  This  puts  the  expression  of  centrifugal  effects  in  the  same 


format  as  the  effects  due  to  internal  strain  energy.  In  other  words,  both  become  structural 


stiffness  terms.  The  internal  strain  energy  components  are  represented  as  conventional 


infinitesimal  small-displacement  elastic  stiffness.  The  components  due  to  centrifugal 


potential  are  called  differential  stiffness.  The  total  of  the  elastic  and  differential  stiffnesses 


is  called  the  tangential  stiffness.  The  proof  that  this  method  accurately  represents  the 


stiffening  effect  of  centrifugal  force  is  in  the  next  chapter  in  the  examination  of  the 


Southwell  effect. 


The  Coriolis  terms,  being  a  function  of  velocity,  cannot  be  represented  as  a  potential. 


Instead,  the  kinetic  energy  formulation  must  be  used.  The  centrifugal  and  Coriolis  effects 


therefore  are  formulated  separately  in  this  analysis.  This  is  permissible  because  the 


Coriolis  acceleration  terms  and  the  centripetal  acceleration  terms  (those  that  generate 


centrifugal  forces)  add  linearly  (Greenwood,  1965).  Therefore,  each  effect  will  be 


considered  separately  and  then  combined  to  produce  the  total  desired  effect. 


For  the  generalized  beam  element,  each  node  has  three  translational  degrees  of 


freedom  and  three  rotational  degrees  of  freedom.  For  the  nodal  displacement  vector  {q} 


given  by: 


{q)T  =  (xx  yi  zi  *X1  9yx  9V  x2  y2  z2  ex2  ey2  ez2/ 


(3.35) 


where  8  is  the  rotation  degree  of  freedom  for  the  indicated  direction,  the  element  matrices 


will  have  dimension  12.  Therefore,  the  element  Coriolis  matrix  will  be  a  12x12  with  all 


forms  equal  to  zero  except  for  the  following  nonzero  components  by-. 


b12  =  2mf) 


b2  j  =  -2mQ 


b78  =  2mQ 


b8  7  =  -  2mfl 
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Once  again,  it  is  noted  that  the  element  Coriolis  matrix  is  skew  symmetric. 

It  is  important  at  this  time  to  note  the  effect  that  the  form  of  each  matrix;  [M],  [B], 
and  [K],  on  system  stability.  The  kinetic  energy  term  T2  by  definition  is  a  positive  definite 
or  positive  semidefinite  function  of  the  generalized  velocities.  Therefore,  the  mass  matrix 
must  always  be  positive  definite  or  positive  semidefinite  (Meirovitch,  1980). 

Generalizations  about  the  matrix  [K]  revolve  around  the  nature  of  the  eigensolutions.  A 
positive  definite  stiffness  matrix  will  result  in  real,  positive  eigenvalues  for  the  system 
solution.  A  positive  semidefinite  stiffness  matrix  will  also  give  some  zero  eigenvalues  with 
the  rest  of  the  eigenvalues  being  real  and  positive.  Real,  nonnegative  eigenvalues  imply  a 
stable  system;  therefore  a  positive  definite  or  positive  semidefinite  stiffness  matrix  will 
result  in  a  stable  system.  The  semidefinite  case  will  admit  the  existence  of  rigid  body 
modes.  If  the  stiffness  matrix  is  negative  definite,  negative  semidefinite,  or  indefinite, 
negative  eigenvalues  will  result.  This  will  generate  a  divergent  eigensolution  and  an 
unstable  system. 

The  Coriolis  matrix  [B]  results  from  the  kinetic  energy  term  Tx.  By  the  definition  of 
Tj,  the  Coriolis  matrix  is  always  skew  symmetric.  However,  Coriolis  forces,  also  known 
as  gyroscopic  forces,  do  no  work  on  the  system.  Therefore,  a  system  with  Coriolis  effects 
not  considered  that  is  stable  and  energy  conserving  will  remain  so  with  the  effects 
included.  If,  however,  a  system  has  damping,  the  structural  matrix  resulting  from  these 
effects  will  be  positive  definite  or  positive  semidefinite.  These  forces  serve  to  dissipate 
energy  from  the  system  and  therefore  the  system  will  no  longer  be  energy  conserving. 
However,  it  will  still  be  stable. 

Summary 

We  have  seen  that  the  problem  of  the  free  vibration  of  a  rotating  system  that 
encorporates  both  centrifugal  and  Coriolis  effects  will  involve  complex  modal  analysis.  The 
eigenvalues  will  be  pure  imaginary;  therefore  the  frequencies  will  be  real.  However,  the 
eigenvectors  wifi  be  complex  with  complex  amplitude  ratios.  This  implies  that  for  a  given 


I, 


eigenvector,  the  eigenvector  components  will  oscillate  out  of  phase  but  at  the  same 
frequency,  i.e.,  the  motion  will  not  be  synchronous. 

The  element  matrices  can  be  derived  by  examining  the  system  energy.  The  Coriolis 
matrix  results  from  the  expression  of  kinetic  energy  that  incorporates  the  mixed  product  of 
nodal  displacement  and  nodal  velocities.  In  the  expression  for  the  free  vibration  of  the 
rotating  system  subject  to  centrifugal  and  Coriolis  effects,  the  Coriolis  matrix  has  terms 
that  multiply  the  unknown  velocity  components.  Since  the  Coriolis  matrix  is  skew 
symmetric,  unlike  a  symmetric  damping  matrix,  the  Coriolis  matrix  is  energy  conserving. 
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CHAPTER  4 


CENTRIFUGAL  EFFECTS  ON  VIBRATION  OF 
ROTATING  BEAMS  AND  BLADES 

The  effect  of  centrifugal  forces  on  the  behavior  of  a  rotating  structure  may  be 
analytically  modeled  by  examining  in  detail  the  centrifugal  terms  in  the  energy  expression. 
One  approach  is  to  place  these  terms  in  the  kinetic  energy  expression  and  obtain  a 
centrifugal  effects  matrix.  This  has  been  done  by  Subrahmanyam  and  Kaza,  Olhoff  and 
Parbery,  and  others.  Another  approach  is  to  create  a  potential  energy  expression  that 
generates  the  centrifugal  effect  and  use  this  expression  to  obtain  a  new  stiffness  matrix 
called  differential  stiffness.  This  approach  is  taken  by  Hoa  and  by  MacNeal.  The  purpose 
of  this  chapter  is  to  prove  the  validity  of  this  approach.  This  will  be  done  by  comparing 
the  analytical  results  using  differential  stiffness  to  the  empirical  Southwell  approach.  The 
method  will  also  be  compared  to  results  obtained  by  Isakson  and  Eisley  using  the  kinetic 
energy  approach  to  centrifugal  effects. 

General  Problem  Statement 

In  the  total  optimization  approach,  the  forward  problem  encompasses  the  modal 
analysis  of  the  rotating  blade.  This  can  be  thought  of  as  a  free  vibration  problem,  with  a 
centrifugal  force  applied  to  the  blade  acting  as  a  stiffening  effect.  This  stiffening  effect  can 
be  quantified  by  examining  the  nonlinear  components  of  the  structural  stiffness  matrix 
(Appendix  B). 

In  nonlinear  analysis,  [KT]  is  the  tangent  matrix: 

[Kt]  =  [K0]  +  [Kd]  +  [Kl]  (4.1) 

where 
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[Kt]  is  the  tangent  stiffness  matrix, 

[K0]  is  the  small  displacement  stiffness  matrix, 

[Kd]  is  the  differential  stiffness  matrix 
(often  called  geometric  stiffness  or  initial  stress  matrix), 

and 

[K-l]  >s  large  displacement  matrix. 

In  solving  the  rotating  blade  problem,  it  is  assumed  that  the  displacements  are 
small;  therefore,  [KL]  is  neglected.  The  differential  stiffness  matrix,  which  is  retained, 
does  not  explicitly  contain  displacements,  but  is  dependent  on  the  stress  level  (Zienkiewicz, 
1977).  Dropping  the  large  deflection  effects: 

[KTqt]  =  [Kg]  [Kq]  (4.2) 

This  total  stiffness  matrix  incorporates  the  load  dependent  terms  and  can  now  be  used  in 
the  free  vibration  problem  to  find  frequencies  and  mode  shapes. 

The  matrix  [K0]  is  given  by: 

[Kq]  =  /v[B]T[D][B0]dV  (4.3) 

where 

[B0]  is  the  strain  shape  function  matrix, 

[D]  is  the  elastic  strain-displacement  matrix, 
and  dV  is  the  incremental  volume. 

Furthermore,  it  can  be  shown  (Appendix  B)  that: 

[Kd]{Au}  =  Jv  ( ABl]T{a)d  V  *4.4) 

where 

[ABl]  is  the  incremental  change  of  the  tangent  strain  matrix  due  to  a  small 
increment  in  the  displacements  and  forces 
and  {a}  is  the  stress  vector. 

Let  us  consider  the  vibration  of  a  simple  model  for  a  first-stage  turbine  blade  on  a 


circular  hub  in  a  typical  aircraft  engine,  such  as  the  General  Electric  CF6  (Figure  2).  The 
blade  outer  radius  is  about  787.4  mm,  the  inner  radius  is  393.7  mm,  the  blade  width  is 
556.8  mm,  and  the  blade  thickness  is  87.38  mm.  The  structure  is  made  of  Inconel  718, 
with  E  (Young’s  modulus)  equal  to  202696  MPa,  p  (density)  equal  to  8.2121E-9  Mg/mm3, 
and  v  (Poisson’s  ratio)  equal  to  0.29.  The  blade  is  rotated  at  a  constant  speed,  Cl. 


Several  rotational  speeds  were  applied,  and  the  first  flexural  mode  was  examined. 
Aj  -  A° 

A  log-log  plot  of - vs.  Q  was  made,  where  A,  and  A?  are  the  fundamental  (first) 

1° 

-*1 

eigenvalues  of  the  rotating  and  nonrotating  structures,  respectively  (Figure  3).  The  valid 

portion  of  the  curve  is  for  Cl  above  40.0  sec- 1.  For  Cl  less  than  40.0  sec-  ^  numerical 

\  _  \0 
Ai  Ai 

errors  give  erroneouslv  high  values  for - .  A  conventional  Cartesian  plot  (Figure  4) 

A° 

A1 

helps  put  these  errors  in  perspective. 


Southwell  Coefficient 

The  Rayleigh-Southwell  approximation  is  an  analytical  method  of  calculating  and 
simply  representing  the  effect  of  centrifugal  stiffening  on  modal  frequency  (Isakson  and 
Eisley,  1960).  Using  this  method,  the  natural  frequency  under  rotation  may  be  expressed 
as: 

4„  =  “»Nn  +  Kn^2  '4.5) 

where 

u/Rn  =  natural  frequency  under  rotation  for  the  nth  mode. 

=  natural  frequency  with  no  rotation  for  the  nth  mode,  and 

Kf,  =  the  Southwell  coefficient  for  the  nth  mode. 

For  plate-like  bodies,  an  analysis  of  the  syst.  n  kinetic  energy  (Isakson  and  Eisley, 
1960)  demonstrates  that  the  Southwell  coefficient  can  be  shown  to  be: 


Figure  2.  Rotating  Blade  Model 

where 


Ti  =  /^xdx 

and  are  displacements  in  the  y-  and  z-directions  describing  the  shape  of  the  nth 

natural  mode  of  the  nonrotating  blade,  and  R  =  rotor  radius. 

Isakson  and  Eisley  in  1960  and  in  1964  computed  the  Southwell  coefficient  for 


various  blade  configurations,  using  analytical  and  experimental  techniques.  However,  a 
careful  interpretation  of  Figure  3  will  provide  a  more  straightforward  method  or 
determining  the  Southwell  coefficient  for  a  particular  mode. 


Figure  3.  Effect  of  Rotation  on  ! 

System  Fundamental  Eigenvalue  ’ 

(Ax  -  A°)/A° 


The  eigenvalue  A  is  related  to  the  frequency  w  by: 
A  =  u.'2 

Therefore,  equation  (4.5)  may  be  rewritten  as: 

An  =  +  Knf?2 


14.7) 


14.81 


where  An  =  eigenvalue  of  the  rotating  structure. 

Rewriting  and  applying  to  the  fundamental  mode  (n=  1)  results  sn: 
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Figure  4.  Cartesian  Plot  of  Rotational  Effect 
on  Nondimensionalized  frequency,  (Aj  -  A°)/A° 
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A  —  A0 

Aj  A  j 

Equation  (4.9)  shows  that  the  plot  of - vs.  Q  should  pass  through  the  origin. 

This  is  demonstrated  in  Figure  4,  with  only  a  small  numerical  error  evident.  Thus  it  is 
seen  that  the  deviation  of  Figure  3  from  a  straight  line  at  low  values  of  Q  is  indeed  caused 
by  small  numerical  errors. 


Taking  the  log  of  Equation  (4.9)  gives: 
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In 


Ki 


+  In  n2 


(4.10) 


The  numerical  curve  in  Figure  3  may  be  approximated  by: 
Ai  -  Afl 


In 


=  ln/3  +  alnQ 


(4.11) 


where  a  and  0  are  unknown  parameters.  The  slope  a  of  the  log-log  curve  is  given  in 
Equation  (4.11)  and  it  is  the  power  to  which  Q  is  raised.  Comparing  Equation  (4.11)  to 
Equation  (4.10)  one  sees  the  theoretical  value  of  a  equals  2.  The  numerical  curve  in 
Figure  3  has  a  slope  of  1.91.  This  indicates  an  error  in  the  finite  element  approximation  of 
4.5%  for  the  experiment. 

Figure  5  illustrates  that  a  may  be  approximated  by  2.  Therefore,  Equation  (4.10) 
may  be  rewritten: 


In 


=  ln/3  +  InQ2 


(4.12) 


Comparing  Equation  (4.12)  to  Equation  (4.10),  one  obtains: 

Ki 

0  —  —  (4.13) 

A0 

Al 

Therefore,  Equation  (4.10)  may  be  rewritten: 

h  ~  A? 

-  =  /3ft2  (4.14) 

A0 

Al 

Ax  -  A® 

This  implies  that  0  is  the  slope  of  a  plot  of - vs  ft2  (Figure  5).  For  this  example.  3 

A0 

Al 

equals  8.0761E-5. 

is  given  from  Equation  (4.13): 

Kx  =  /3A°  (4.15) 


For  this  problem,  Kj  =  46.45. 
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z-displacements  (no  y-rotations). 


Figure  6.  Rotating  Beam  with  End  Mass 


The  extensional  (x-direction)  displacements  are  needed  to  generate  internal  loads 


(were  these  displacements  not  permitted,  the  load  P  would  be  totally  absorbed  by  the 


constraint).  There  is  no  differential  stiffness  in  the  extensional  direction  (Appendix  B). 


The  extensional  components  of  the  stiffness  matrix  are  not  affected  and  the  eigenvalue  in 


this  direction  remains  unchanged.  Therefore,  for  the  purpose  of  this  demonstration,  we 


will  concern  ourselves  only  with  masses,  stiffnesses,  and  eigenvalues  involving  only 


transverse  (z-direction)  displacements.  This  can  be  done  because  the  stiffness  and  mass 


terms  for  the  x-  and  z-directions  are  completely  uncoupled. 


We  now  reduce  the  small  displacement  stiffness  matrix  [K0]  to  one  term: 


[Kq]  = 


(4.17) 


The  differential  stiffness  matrix  [KD]  has  nonzero  components  of  only  one  term  and  is: 


(4. 


where  P  is  the  load  applied  to  the  beam.  In  this  case,  P  is  the  centrifugal  force  and  is 
given  by: 

P  =  4ir2n2m L  (4. 


Substituting  Equation  (4.19)  into  Equation  (4.16)  and  using  Equation  (4.2),  one 
obtains: 

12EI  24ir2n2m 

[K-TOtJ  =  -  4-  — — —  (4. 
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The  general  form  of  the  eigenvalue  problem  in  finite  element  analysis  is: 

[Krpo'rH'f1}]  =  Aj[M]{0}j  (4. 

where  [M]  is  mass  matrix,  {<p}l  the  eigenvector,  and  A;  the  eigenvalue. 

The  eigenvalues  Aj  are  found  from  solving  the  characteristic  equation: 

I  I^-TOT^  ~  \[M]  |  =  0  (4. 

This  particular  problem  considers  only  one  degree  of  freedom;  therefore,  only  one 
eigenvalue  Ax  exists.  It  is  found  from  Equation  (4.20)  to  be: 

12EI  24t2Q2 

Ax  =  -  +  -  (4. 
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The  eigenvalue  in  the  z-direction  for  the  nonrotating  problem  A°  can  be  found  by 
using  the  stiffness  matrix  given  in  Equation  (4.17): 

12EI 

Aj  =  -  (4. 

mL3 


(Note  that  the  above  result  could  also  be  derived  from  Equation  (4.23)  with  Q=0.) 
Equations  (4.23)  and  (4.24)  can  be  manipulated  to  obtain: 


Comparison  of  Equation  (4.25)  to  Equation  (4.11)  obtains: 


5 

0  =  - 

El 


Therefore,  Equation  (4.13)  gives: 


(4.26)  ■ 

t 


Kj  =  ^r2  (4.27)  ; 

~  47.3741  f 

) 

{ 

This  is  the  theoretical  value  of  the  Southwell  coefficient  for  the  mass  with  one  degree  of  j 

freedom  at  the  end  of  a  massless  rotating  beam. 

Let  us  now  try  a  numerical  experiment.  Let  E  equal  73.77E3  MPa,  L  equal  250 

mm,  I  equal  3.2552E4  mm4,  A  (cross  section  area  of  the  beam)  equal  625  mm2,  and  m 

\  —  \0 
Ai  Ai 

equal  5.0E-4  Mg.  Figure  7  is  a  log-log  plot  of - vs.  Cl  and  Figure  8  is  a  conventional  ■ 

i° 

1  Al_A?  02  ! 
Cartesian  plot  of  the  same  variables.  Figure  9  is  a  plot  of - vs  Cr.  The  finite  j 

i° 

Ai 

element  data  agree  exactly  with  the  theoretical  results  given  by  Equation  (4.25).  Figure  7  ! 

shows  that  a  =  2.000.  0,  derived  from  Figure  9,  =  1.2844E-5.  Equation  (4.15)  gives  Kx 
=  47.374.  This  agrees  exactly  with  the  theoretical  result  of  Equation  (4.27).  Therefore, 

t 

for  this  problem,  there  is  no  error  between  the  finite  element  approximation  and  the  J 

theoretical  result.  [ 

Thus  we  see  that  finite  elements  and  graphical  analysis  can  provide  a  simpler 
method  for  calculating  the  Southwell  Coefficient.  This  technique  is  more  general  than  that  j 

of  Isakson  and  Eisley  in  that  they  were  restricted  to  blade  shapes  that  could  provide  a 
closed  form  solution.  The  finite  element  method  can  provide  Southwell  type  of  information 
for  a  rotating  blade  of  any  geometry.  j 

It  is  also  seen  that  the  potential  energy  approach  that  leads  to  the  development  of  < 

1 

¥ 

differential  stiffness  can  be  used  to  accurately  determine  the  stiffening  effect  of  rotation  on  t 

the  eigenvalues  of  the  structure.  Results  obtained  from  this  method  compare  favorably  j 

with  those  obtained  from  a  theoretical  examination  of  the  change  in  the  characteristic  | 

I 
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CHAPTER  5 


PERTURBATION  METHODS  IN  OPTIMAL  DESIGN 

This  chapter  will  discuss  the  theory  and  applications  of  the  perturbation  method  as 
it  relates  to  the  problem  of  redesign  in  the  presence  of  rotational  effects.  The  theory  of 
nonlinear  inverse  perturbation  will  be  presented.  The  linear  perturbation  equations  as 
derived  by  Stetson  will  be  modified  to  include  centrifugal  effects.  Then,  new  equations  for 
nonlinear  inverse  perturbation  will  be  derived  from  first  principles  that  include  complex 
components.  These  equations  will  be  applicable  for  Coriolis  effects.  An  overview  of  the 
optimization  techniques  used  will  also  be  presented. 

Definition  of  Basic  Terms  in  Optimization 

The  baseline  system  is  the  original  system  configuration.  The  starting  values  of  the 
system  parameters  and  system  performance  are  first  determined  for  the  baseline  system. 

The  objective  system  is  the  system  that  exhibits  the  desired  system  performance. 

This  system  is  the  goal  of  the  redesign  process. 

The  system  parameters  are  those  independent  factors  that  each  influence  the  system 
performance.  One  such  parameter,  and  the  one  that  will  be  the  primary  subject  of  this 
analysis,  is  element  thickness.  This  parameter  influences  both  mass  and  stiffness.  One 
goal  of  the  first  part  of  the  analysis  is  to  determine  how  changes  in  thickness  actually  alter 
the  mass  and  stiffness  matrices. 

The  system  performance  is,  simply,  the  behavior  of  the  system.  System  performance 
can  be  measured  in  terms  of  statics,  such  as  stress  and  displacement,  or  in  terms  of 
dynamics,  such  as  eigenvalues  and  eigenvectors. 

Optimal  redesign  is  the  process  by  which  the  desired  change  in  performance  is 


obtained  by  altering  the  parameters  in  the  most  advantageous  manner  so  as  to  minimize 
(or  maximize)  the  system  objective  function.  This  objective  function  is  a  measure  of  a 
desired  quantity.  Two  of  the  more  common  system  objective  functions  are  minimimum 
weight  and  minimum  change  from  the  baseline  system. 


Perturbation  Methods 

The  first  task  is  to  determine  the  effect  of  the  system  parameters  on  system 
performance.  This  involves  analyzing  the  structural  matrices  to  see  what  the  effect  of  a 
change  in  one  basic  parameter,  thickness  for  example,  would  have  on  each  matrix.  This 
having  been  done,  the  effect  of  changes  to  the  baseline  system  can  be  tracked  throught  the 
redesign  process  as  one  works  toward  the  objective  sj^stem. 

We  first  wish  to  determine  the  influence  of  element  thickness  upon  the  differential 
stiffness  matrix  of  beam  and  plate  elements.  This  matrix  may  be  defined  by  (Cook, 

1974): 

t  „  [s]  [0]  [0] 

[Kfj]  =  /  [N']t  [0]  [s]  CO]  [N']dV  (5.1) 

Jv  [0]  [0]  [s] 

where 

[Kq]  =  elemental  differential  stiffness  matrix, 

[N'3  —  derivative  of  the  shape  function  matrix,  and 
[s]  =  matrix  of  applied  stresses 

axQ  rxy0  rxz0 
[s3  =  Txy0  ffy0  Tyz0 
Txz0  ryz0  °z0 

Figure  2  from  Chapter  4  illustrates  a  typical  rotating  structure.  Let  ff  be  the 
rotational  frequency  of  the  structure  about  the  z-axis  in  hz.  Let  us  now  suppose  that  there 
exists  an  element  of  the  structure  with  the  centroid  at  a  radial  distance  r0  from  the  origin, 
length  b,  mass  m,  thickness  t,  cordwise  distance  a,  and  density  p.  The  centrifugal  force  F 


applied  to  this  structure  is  given  by: 

F  =  4x2n2r0m  (5.2) 

where  m  may  be  expressed  by: 

m  =  pabt  (5.3) 

Therefore, 

F  =  4x2ft2r0pabt  (5.4) 

This  illustrates  that  the  centrifugal  force  is  linearlj7  dependent  on  thickness.  The 
differential  stiffness  matrix  is  linearly  dependent  on  the  centrifugal  force,  though  the 
constant  of  proportionality  may  be  dependent  on  geometry.  Therefore,  we  can  say: 

[Kp]  oc  t  (5.5) 

This  implies  that  a  change  in  thickness  will  effect  the  differential  stiffness  matrix  linearly. 

One  can  relate  the  baseline  and  objective  systems  through  perturbations  of  the 
baseline  system  quantities.  The  stiffness  and  mass  matrices  are  perturbed  by: 

[kl  =  [k]  +  [Ak]  (5.6) 

and 

[m'j  =  [m]  +  [Am]  (5.7) 

where  [Ak]  and  [Am]  are  the  perturbations  to  the  baseline  stiffness  and  mass  matrices, 
respectively.  These  perturbations  will  cause  perturbations  in  the  dynamic  response.  Let 
the  matrix  of  eigenvectors  be  given  by  [<£]  and  the  diagonal  matrix  of  eigenvectors  be  given 
by  [w2].  The  perturbations  in  the  eigenvalue  and  eigenvector  matrices  are  given  by: 

[u/2]  =  [a-2]  +  [  A(w2)]  (5.8) 


W  1  -  [ <t> ]  +  [  A<j>] 


where  [A(w2)]  and  [A$]  are  the  changes  to  the  baseline  system  eigenvalues  and 


eignevectors,  respectively. 


The  structural  changes  described  in  Equations  (5.6)  and  (5.7)  can  be  decomposed 
into  p  element  change  properties  where  a  group  of  elements  may  be  allowed  to  change. 


*  *  r  t  **  »«»L  >»4rf4Vi  tlfcVI  tU«Vj  4*4,u  Jlj', 


[6k]  =  £[Ake]  (5.10) 

e  =  1 


[Am]  =  £[Ame]  (5.11) 

e  =  1 

where  [Ake]  and  [Ame]  are  changes  to  the  stiffness  and  mass  matrices  of  element  e, 
respectively.  Furthermore,  each  element  change  can  be  expressed  as  a  fractional  change 
q  from  the  baseline  system.  The  change  a  may  represent  a  change  in  element 


thickness.  In  general,  ae  can  be  expressed  by: 


[Ake]  =  [ke]ae 


[Ame]  =  [me]ae 


(5.12) 


(5.13) 


Hoff  (1985)  used  the  linear  relationships  given  in  Equations  (5.12)  and  (5.13)  in  his 
work  on  optimal  redesign;  however,  he  pointed  out  that  nonlinear  relations  may  be 
required  for  certain  applications.  In  the  case  of  plates,  membrane  components  of  the 
stiffness  matrix  and  also  the  differential  stiffness  matrix  vary  linearly  with  the  plate 
thickness.  The  bending  component  of  the  stiffness  matrix,  however,  varies  as  the  cube  of 
the  plate  thickness.  Therefore,  for  plates,  Equation  (5.12)  is  replaced  by: 


[Ake]  [kemembJQe  + 


+  [k«bend](3Qe  +  3q'*  + 


(5.14) 


where  tkememb^  contains  the  membrane  components,  [kedinJ  contains  the  differential 
stiffness  components,  and  [k„,  ,]  contains  the  bending  components  of  (kj. 

ebend  c 

Equation  (5.14)  also  holds  for  beams,  with  the  exception  that  the  element  stiffness 

matrix  containing  only  extensional  properties.  fkeext^  reP*aces  fkemembl-  Therefore,  for 


beams: 


[Ake]  [keextlQe  +  ^kedifpQe 


+  [kebendK3Qe  +  3a*2  +  a*3) 


(5.15) 


The  first  matrix  eigenvalue  and  eigenvector  redesign  method  using  perturbation  was 
developed  by  Stetson  (1975).  This  procedure  involves  the  linearization  of  the  uncoupled 
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energy  equations  of  the  objective  structure: 

[K'l  =  (5.16) 

Additional  work  on  this  method  was  done  by  Stetson  and  Palma  (1976),  Stetson  et 
al.  (1978),  Sandstrom  and  Anderson  (1981),  and  Hoff,  Bernitsas,  and  Anderson  (1985). 

Stetson  developed  the  generalized  form  of  the  perturbation  equation: 

[0']T[Ak]m  -  [<*’]T[Am]  [4>']W2]  =  t^']T[m][^'][w'2]  -  [0]T  [k][0']  (5.17) 

To  solve  the  above  equations,  note  that  Equation  (5.17)  is  composed  of  n2  scalar  equations 
given  by: 

W'}JT[Ak]{^'}I  -  Wi’2OTj[Am]  ~  WjTO'h  (5.18) 

for  ij  =  1,2,... ,n. 

The  changes,  Am  and  Ak,  were  defined  in  Equations  (5.10)  and  (5.11).  The  mode 
shape  changes  [A0]  were  given  by  Stetson  in  his  first  order  perturbation  by: 

[A*]  =  [0][c]T  (5.19) 

where  the  admixture  coefficients  Cy,  i,  j  =  1,2,.. .n  are  small  and  cu  =  0. 

Equations  (5.8)  and  (5.9)  tire  applied  to  Equation  (5.17).  This  results  in: 

[4>f[  AkM  +  2[0]T[Ak][A^>]  +  [  A*]t[  Ak][  A<tf 

H[0]T[Am][0]  +  2[0]T[Am][A0]  +  [A0]T[Am][A0])  ([w2]  +  (A(u.'2)]) 

(5.20) 

=  ([0]T[m][d>]  +  2[0]T[m][  A0]  +  [  A<«T[m][ A*])  ([u,2]  +  [  A(w2)]) 

-  (MT[k]M  +  2(0]T[k][  A0]  +  [  A0]T[k][  A<«) 

Stetson  then  obtains  the  first  order  perturbation  equations  by  neglecting  terms  of  order  A2, 
A3,  and  A4.  This  yields  (in  scalar  form): 

{y>}J[ Ak){v}j  “  {vjJtAmKw}^,2  =  M^A^2)  for  i  =  j  (5.2 1 ) 

=  MjCy(u>j2  -  Wj2)  for  i  *  j 

The  solution  of  the  above  first  order  equations  require  the  specification  of  the 
frequency  changes  A(wj2)  and  the  mode  shape  changes  Awk]  in  terms  of  the  admixture 
coefficients  c,i  where  Apkl  is  the  change  in  the  kth  degree  of  freedom  of  the  ith  mode. 

In  order  to  eliminate  the  admixture  coefficients,  the  following  algebraic 
manipulations  are  performed  (Sandstrom  and  Anderson,  1982)  for  the  case  of  non- 


repeated  eigenvalues.  The  change  to  the  kth  degree  of  freedom  of  the  ith  mode  in  terms  of 
the  admixture  coefficients  is: 


^ki  =  cil^kl  +  ci2^k2  +  -  +  ci«An 

n  (5.22) 

=  Y 

J=1  ,  j*i 

Rearranging  Equation  (5.21)  for  i*j  obtains  the  following  expression  for  the  admixture 


coefficients: 


1 

c  - - Ak]{i/>}j  -  wi2{V’}J[Am]{V<}i)  (5.23) 

M^2  -  Wj2) 

Applying  Equation  (5.23)  to  Equation  (5.22)  results  in  the  following  expression  that  relates 
the  physical  mode  shape  changes  directly  to  the  structural  changes: 


n 


ki  =  Y 

j=i 


_Mj(uj2  -  Wj2) 


(5.24) 


for  i*j. 

Using  the  relationship  for  the  change  in  the  element  mass  matrix,  Equation  (5.14), 
and  the  nonlinear  relationship  for  the  change  in  the  element  stiffness  matrix,  Equation 
(5.15),  results  in  the  following  expression,  nonlinear  in  the  element  change  property  ae, 
that  describes  the  change  in  the  natural  frequency  to  the  ith  mode: 

+  nWMft 

i  e=  1 

+  {^tkebend]Wi(3ae  +  ^  +  f5’25) 

-  ^,2{w}J[me]  {v},ae 

Also,  the  nonlinear  expression  in  terms  of  ae  for  the  physical  mode  shape  change  for 


the  kth  degree  of  freedom  becomes: 
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^ki=  I  E  ~  ~T  {{^([kememb]  +  rkediftJW-Qe 

e  =  i  i,J =  1  —  Wj  )  L 

+  {^}JlKbend^h(3ae  +  3ae2  +  ae3) 


(5.26) 


-  1  (for  j 


In  applying  the  method  described  above  in  finite  element  analysis,  practical 
considerations  make  it  necessary  to  divide  the  quadrilateral  elements  in  the  finite  element 
model  into  two  elements:  one  with  only  membrane  stiffness  and  one  with  only  out-of-  plane 
(bending)  stiffness.  These  elements  are  then  superimposed.  This  permits  multiplication  of 
the  stiffness  terms  representing  membrane  properties  by  a  linear  element  change  factor 
while  the  stiffness  terms  containing  the  bending  properties  can  be  altered  by  a  nonlinear 
change  factor. 


Derivation  of  Complex  Perturbation  Equations 
Let  us  now  determine  the  changes  in  the  system  eigenvalues  and  mode  shapes  for  a 
problem  exhibiting  Coriolis  behavior.  As  discussed  previously,  the  eigenvalues  for  this 
system  will  be  real;  however,  the  mode  shapes  will  have  complex  components. 


The  basic  equation  of  motion  for  the  descretized  system  is: 

[MK&i  +  [BMj  +  [K]{^}j  =  {0} 


(5.27) 


where  [M],  [B],  and  [K]  are  the  system  mass,  Coriolis,  and  stiffness  matrices  respectively, 
and  {<b}j  is  the  displacement  vector  for  the  ith  mode. 

Let  us  assume  a  solution  for  the  system  described  in  Equation  (5.27)  of  the  form: 


Wh  = 


(5.28) 


where  {^},  is  the  eigenvector  for  the  ith  mode  and  is  the  eigenvalue  for  the  ith  mode. 


Applying  Equation  (5.28)  to  Equation  (5.27)  obtains: 

-^[MlMie^+i^fBlW^1  +  [k]{*}iei“it  =  {0} 


(5.29) 


The  term  e"'1  will  be  eliminated.  In  addition,  the  terms  will  be  premultiplied  by  Wf, 
to  obtain  the  following  scalar  equation: 
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+  Wi+WflKlW-0  (5.30) 

Define  generalized  mass  for  mode  i,  M;,  to  be: 

Mi  =  {^[M]Mi  (5.31) 

Similarly,  generalized  damping  or  generalized  Coriolis,  Bj,  is  defined  by: 

Bi=W?[B]Wi  (5.32) 

Final!}',  generalized  stiffness  may  be  defined  as: 

Ki={^[K3Wi  (5.33) 


Stetson  (1975)  applied  the  above  definitions  for  generalized  mass  and  generalized 
stiffness  in  obtaining  his  perturbation  equations.  Since  the  following  analysis  will  employ 
both  real  and  imaginary  components,  the  full  form  of  the  terms  in  Equation  (5.30)  will  be 
retained  for  now. 

Perturbations  of  the  System  Including  Complex  Effects 

When  the  original  system  is  modified  in  the  optimization  process,  it  can  be  said  to  be 
perturbed.  The  perturbed  system  must  also  obey  the  equations  of  equilibrium.  Let  the 
perturbed  system  be  distinguished  from  the  original  by  primes.  Equation  (5.30)  can  be 
rewritten  for  the  perturbed  system  by: 

-W]T[M’]m[«'2]  +  i[1/«,]T[B,]W][w']  +[^]T[K'][^']  =  0  (5.34) 

where  is  the  matrix  of  perturbed  eigenvectors,  and  [V]  is  the  diagonal  matrix  of 
perturbed  eigenfrequencies. 

The  perturbed  system  can  be  related  to  the  original,  unprimed  system  by: 


[K'l  =  [K]  +  [AK] 

(5.35) 

[M']  =  [M]  +  [AM] 

(5.36) 

[B'l  =  [B]  +  [AB] 

(5.37) 

[tl/]  =  [w]  +  [Aw] 

(5.38) 

[w’2]  =  [w2]  +  [Aw2] 

(5.39) 

W]  =  [v]  +  [Ai/j] 

(5.40) 

The  above  perturbation  equations  are  similar  to  Equations  (5.7)  through  (5.10). 
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Equations  (5.11)  and  (5.12)  define  the  global  changes  in  mass  and  stiffness  in  terms  of 
element  change  properties  and  Equations  (5.13)  through  (5.15)  define  the  element  changes. 
One  last  change  that  needs  to  be  defined  is  that  for  the  Coriolis  matrix.  This  is  done  by: 


[AB]  =  £[Abe] 

e=  1 


[AbJ  =  [bja 


(5.41) 


(5.42) 


The  last  equation  is  justified  because,  as  we  have  seen,  the  Coriolis  matrix  is  linearly 
dependent  on  element  mass.  The  element  mass  is  itself  linearly  dependent  on  changes  in 
thickness;  the  element  change  parameter. 

Equations  (5.35)  through  (5.42)  are  applied  to  (5.34).  Terms  in  A  of  order  2  or 
higher  are  eliminated  as  are  the  baseline  equilibrium  terms,  which  essentially  are  order 
zero  in  A.  In  a  manner  similar  to  that  used  in  section  4.1,  the  equations  are  expanded  for 
all  modes  and  an  expression  for  the  change  in  the  eigenvalue  in  terms  of  admixture 
coefficients  is  obtained.  This  relates  the  change  in  the  eigenvectors  for  all  modes  to  the 
changes  in  element  properties  during  the  redesign  process.  This  results  in  the  following 
expression  for  the  change  in  natural  frequency  in  scalar  form: 


{^/[AKJMi  +iwi{V’}/[AB]{V'}i  -u>f{ifi}j[  AM]{V>}j 


=  M?[MMjA (uf)  (i=j)  (5.43) 

=  (Wj[M]{V»}j(w^-wf)  —  i{V'}j'[B]{-0}i(  Atoj  —  Atoi))cij  (i*j) 

Applying  the  definitions  for  the  changes  in  the  structural  matrices  to  the  above 
equation  and  setting  i=j,  the  following  expression  is  obtained  for  the  frequency  change  for 
the  ith  mode  due  to  application  of  the  element  changes  ae: 


,,2'i  -U./,\T 


(5.43) 


{*}? [MMj  A (wf )  -  iW? (BlWj  Awj 

=  £  ({^[k;M,  +iu,i{v}?[beMi  -wfW^nrgWiK 

e=  1 


(5.44) 


where  [k~]  is  the  approximation  of  the  cubic  expansion  of  the  element  stiffness  matrix  such 


that: 


gsnaww^yumaaMMAH  »  m  mthmhh  uwh  ww  lvw  m  w  w  iwuuwjuwwwmwwwwiiwi  uwinjunumimiinB 
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^  =  fkememb]  +  [kJ  +  3[kebend]  (5'45) 

However,  it  must  be  remembered  that  the  eigenvector  is  itself  a  complex  quantity, 

with  both  amplitude  and  phase  components.  Let  us  express  the  complex  eigenvector  for 

the  ith  mode,  {ip}^  as: 

Mi  =  Mi  +  i{£h  (5.46) 

where  {0};  denotes  the  real  part  of  the  ith  eigenvector  and  {£}>  denotes  the  imaginary  part. 
Therefore,  Equation  (5.44)  becomes  two  equations;  one  that  equates  real  parts  and  another 
one  that  equates  imaginary  parts.  These  are,  respectively, 

(MftMMi-tttfCMMi )  A(W?)  +»f[BMGlA(«i) 

p 


=  E{{^fke]Mi  -{OjlKMi  -  -{$[ }a( 


(5.47) 


and 

2{^nM]{^A(u.f)  +{^[B]{e}i)AWi 

p 


=  Ei{2^}i'[k'eMi  +  fwrtbeK^i+WrtbJWijwi  -2M?[meKtfWf  }«, 


(5.48) 


Note  that  for  the  case  of  no  Coriolis  terms  and  purely  real  eigenvectors,  Equation 
(5.47)  degenerates  to  Equation  (5.25)  and  Equation  (5.48)  becomes  identically  satisfied. 

To  determine  the  eigenvector  change  in  the  Coriolis  system,  the  admixture 
coefficients  c^  previously  defined  in  Equation  (5.45)  can  be  used.  From  these  admixture 
coefficients,  the  eigenvectors  of  the  perturbed  system  can  be  obtained: 


P  n 

A^ki  =  E  E 

e=l  j=l 


^kj 


LM^MIMM,2-^2)  -  i{tf}ftB]{*}i(«i-«j) 


(5.49) 


+  iui{i>}J[be]{ii}lae  -  [me]{V'}iae  | 


>(  for  j*i) 


Sequential  Unconstrained  Minimization  Techniques 
The  generalized  form  of  the  optimization  problem  subject  to  both  equality  and 
inequality  constraints  may  be  stated  by: 

min  F({x})  (5.50) 

subject  to 

gj({x})£0  j=  1,2,... ,m  (5.51) 

hk({x})  =  0  k=l,2,...,l  (5.52) 

In  the  external  penalty  function  method,  a  penalty  is  associated  with  the  violation  of 
a  constraint.  The  penalties  are  applied  outside  the  region  of  the  feasible  domain.  All 
intermediate  solutions  lie  in  the  infeasible  region.  The  solution  is  obtained  by  convergence 
from  the  outside.  The  solution  may  start  from  an  infeasible  point;  therefore,  an  initial 
feasible  point  is  not  required.  However,  when  the  optimum  solution  is  achieved,  this 
solution  also  is  not  in  the  feasible  region.  Bellagamba  and  Yang  (1981)  applied  the 
exterior  penalty  function  approach  to  the  problem  of  truss  optimization. 

The  internal  penalty  function  method,  however,  always  keeps  the  solutions  inside 
the  feasible  domain.  The  solution  procedure  can  be  stopped  at  any  time  and  a  permissible 
optimized  result  will  be  obtained.  If  one  wishes  to  start  the  design  process  from  inside  the 
feasible  region,  one  must  obviously  start  with  a  feasible  solution.  This  may  not  always  be 
possible. 

The  Augmented  Lagrange  Multiplier  (ALM)  method  produces  a  way  to  reduce  the 
dependency  of  the  algorithm  on  the  choice  of  penalty  parameters  and  the  way  they  are 
updated.  This  is  accomplished  by  combining  the  use  of  Lagrange  multipliers  with  penalty 
functions.  Using  only  Lagrange  multipliers  gives  an  optimum  that  is  a  stationary  point 
rather  than  a  true  minimum  of  the  Lagrangian.  Using  only  stationary  functions  gives  a 
minimum  that  leaves  open  the  possibility  of  an  ill-conditioned  solution  that  is  not  feasible. 
Using  both  creates  an  unconstrained  problem  that  obtains  a  feasible,  well-conditioned 
solution  that  is  a  true  minimum. 


The  minimization  problem  may  now  be  stated  in  terms  of  the  following  augmented 


Lagrangian: 


F({x},{A},r  )  =  F({x})  +  +  rtf] 

P  M  P  (5.53) 

+  E{^k  +  mhk(W)  +  rp[hk({x})]2} 
k=l  *■  > 


where 


(5.54) 


The  update  formulas  for  the  Lagrange  multipliers  are: 


1 

f 

-Af 

) 

AP+1  =  AP  +  2rp  ( 

j^max 

gj({x}), - 

L  2rP- 

l 

1 

AE  +  m  =  AE  +  m  +  2rphk{x}  {5.56] 

for  j  =  l?m  and  k  =  1,1. 

Vanderplaats  (1984)  states  that  the  ALM  method  has  the  following  advantages: 

1.  The  method  is  relatively  insensitive  to  the  value  of  rp.  It  is  not  necessary  to 
increase  rp  to  oc. 

2.  Precise  gj({x})  and  hk({x})  are  possible. 

3.  Acceleration  is  achieved  by  updating  the  Langrange  multipliers. 

4.  The  starting  point  may  be  either  feasible  or  infeasble. 

5.  At  the  optimum,  the  value  of  A*  *0  will  automatically  identify  the  active 
constraint  set. 


Vj  =  max 


gj{x}, - 

2r„ 


Vanderplaats  (1982)  provides  a  historical  overview  of  optimization  methods. 
Vanderplaats  (1983)  introduced  a  computer  implementation  of  optimization  methods,  ADS 


CHAPTER  6 


OPTIMAL  REDESIGN  METHOD 

In  this  chapter,  the  predictor-corrector  method  for  optimal  redesign  will  be 
developed.  This  development  will  take  place  in  the  context  of  simple  example  problems. 

In  later  chapters  the  method  will  be  expanded  to  large  problems  and  to  the  inclusion  of 
Coriolis  effects  with  complex  eigenvectors. 

Predictor-Corrector  Method 

Consider  the  first  example  for  optimal  structural  redesign,  the  rotating  beam,  shown 
in  Figure  10.  Only  vibration  in  the  x-z  plane  will  be  considered;  therefore,  the  permitted 
vibratory  degrees  of  freedom  will  be  x-  and  z-axis  displacements  and  y-axis 
rotations.  Furthermore,  the  left  end  of  the  beam  will  clamped,  and  the  right  end  will  be 
“guided”  so  as  to  permit  only  x-  and  z-axis  displacements  and  rotations.  This  guided 
bending  condition  is  chosen  because  the  problem  then  can  also  be  readily  solved  by 
theoretical  means.  The  theoretical  solution  will  be  used  for  comparisons. 

The  first  design  change  seeks  a  10%  increase  in  the  fundamental  flexural  model 
frequency.  In  the  predictor  step,  we  will  assume  that  the  element  change  ae  is  small; 
therefore,  the  quantity  (l  +  ae)3  -  1  may  be  approximated  by  3qs.  This  results  in  a 
3.11%  error,  but  is  done  so  as  to  facilitate  solving  for  a,  which  we  shall  see  will  be  the 
unknown  in  the  inverse  perturbation  scheme.  The  element  change  property  oe  is  given  by; 

Qe  =  -  (6.1) 

‘e 

where  te  is  the  element  thickness,  and  Ate  is  the  change  in  element  thickness.  Thickness 
for  the  beam  in  Figure  10  is  z-direction  in  the  cross  section.  The  scalar  equation  that  gives 


Figure  10.  Rotating  Beam 


the  relationship  between  ae  and  the  desired  change  in  natural  frequency  for  the  ith  mode 
Awj,  neglecting  Coriolis  effects  but  including  centrifugal  effects,  is  given  from  Equation 
(5.25): 


^2>  -  w.  E 


'  e=  1 


-  {ip} jQe 


(6.2) 


where  the  approximation  of  the  expansion  of  (l+ae)3  as  3ae  is  applied. 

The  above  equation  is  called  the  predictor.  It  relates  the  change  in  the  element 
change  properties  to  a  prescribed  change  in  the  desired  eigenfrequency.  In  this  way  the 
equation  predicts  what  the  system  configuration  should  be  for  a  given  amount  of  frequency 
change.  For  the  complex  case  involving  Coriolis  effects,  Equations  (5.47)  and  (5.48)  serve 
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as  the  predictor  equations.  Note  that  in  that  case  there  sire  two  predictor  equations. 


For  the  single  element  case,  Equation  (6.2)  can  be  solved  as  one  equation  with  one 


unknown.  When  the  bar  is  divided  up  into  two  or  more  elements,  the  excess  unknown 


element  change  properties  can  be  found  by  optimizing  some  function,  such  as  minimum 


weight  or  minimum  structural  change.  Equation  (6.2)  then  becomes  an  equality  constraint 


in  the  optimization  scheme.  In  the  Coriolis  problem,  there  are  two  equality  constraints. 


Therefore,  it  is  seen  that  the  problem  becomes  one  of  parametric  optimization,  with 


thickness  as  the  parameter  to  be  optimized. 


Using  the  results  from  the  static  analysis,  one  finds  the  eigenvalues  and 


eigenvectors  for  the  rotating  system  by  finite  element  modal  analysis  using  MSC/ 


NASTRAN  as  modified  with  Direct  Matrix  Abstraction  Programming  (DMAP).  The 


restart  procedure  used  is  described  in  Appendix  A.  The  eigenvectors  are  stored,  printed 


out,  and  written  to  a  FORTRAN  file.  Another  FORTRAN  postprocessor  strips  out  the 


eigenvector  for  the  desired  mode  and  appropriately  partitions  it  for  each  element.  Once 


this  has  been  accomplished,  all  of  the  quantities  necessary  to  formulate  Equation  (6.2) 


have  been  obtained. 


Equation  (6.2)  is  used  as  an  equality  constraint  in  the  Augmented  Lagrange 


Multiplier  (ALM)  procedure  described  in  Chapter  5.  Inequality  constraints  are  also 


formulated.  The  first  inequality  constraint  requires  a  to  be  greater  than  —0.5.  This 


ensures  that  the  element  thickness  will  always  be  positive  during  the  redesign  process  and 


in  no  case  will  an  element  be  reduced  by  more  than  50%.  This  makes  certain  that 


unwanted  secondary  effects,  such  as  static  failure  due  to  the  “applied”  centrifugal  load  will 


not  occur.  The  second  inequality  constraint  forces  the  a  to  be  less  than  1E5.  This 


supplies  an  upper  bound  to  the  search  procedure.  The  function  to  be  minimized  in  the  first 


example  is  the  design  change: 


f(W)  =  I>e)2 
e=  1 


I 
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where  {<>}  is  the  vector  of  element  change  properties  ae. 

Alternatively,  the  function  to  be  minimized  could  be  minimum  weight.  For  a  system 
of  uniform  density,  this  function  is  given  by: 

p 

«{<*})  =  £Aeae  (6.4) 

e=  1 

where  Ae  is  the  element  planform  area  perpendicular  to  the  thickness. 

The  element  change  properties  determined  from  Equations  (6.3)  or  (6.4)  are  used  to 
recompute  the  cross  sectional  area  and  moments  of  inertia  for  each  element.  A  reanalysis 
is  then  accomplished  to  determine  the  eigenvalues  and  eigenvectors.  The  perturbed 
eigenvectors  are  necessary  to  perform  the  next  step  in  the  procedure,  namely  the 
corrector. 

The  corrector  examines  the  potential  energy  imbalance  between  the  system  output 
from  the  predictor  and  the  desired  system  and  corrects  the  imbalance  through  additional 
elemental  changes.  This  enforces  the  natural  frequency  constraint  on  the  zth  mode.  The 
following  equation  is  used  for  the  example  with  no  Coriolis  effects: 


e  =  1 


(6.5) 


where  {-0*}  is  the  perturbed  eigenvector,  w'  is  the  desired  eigenvalue,  and  [K]  and  [M]  are 
the  global  stiffness  and  mass  matrices,  respectively.  Note  that  the  global  stiffness  matrix 
is  the  total  tangential  stiffness  matrix  that  includes  the  differential  components. 

If  Coriolis  effects  are  included,  then  the  following  equation  is  used  for  the  corrector: 


+  Wh  -«,?WfrmeHA)ae 

e=l  (6.6) 

=  [BMn  -  WtftWh 

Notice  that  the  above  equation  can  be  broken  down  into  two  equations:  one  that 
equates  the  real  parts  and  another  that  equates  the  imaginary  parts.  This  results  in  two 
equality  constraints  for  the  corrector  step: 


I.*  |.i  •  _«  JU#  (.#  I.*  f.1  { 


•  *.l  h'  tt*  Ik1  1 


!«y5 


£  {^[k-e]{A  -{Oj'Wli  -  [Wjlm.M},  -{n?tme]{niUf  a, 


=  {wjmwh  -{oftMxniV?  -{n^Kiioo 


£  hwfa-M  +lm?[bemi  +{o?n>emX'i  -2 w^my? 


=  2{*,}ftM]{aiu'?  -mjmWi  +{r}?'tB]{f}iy  -2W‘[KHn; 


Notice  that  when  the  Coriolis  effects  are  absent,  the  above  equations  degenerate  into  those 


used  previously. 


The  perturbed  eigenectors  may  be  obtained  in  one  of  two  ways.  The  first  method. 


mentioned  above,  is  simply  rerunning  the  predictor.  This  yields  the  full,  nonlinearly- 


perturbed  matrix  of  eigenvectors  and  the  desired  mode  can  be  easily  partitioned  out.  The 


second  procedure  involves  the  application  of  Equation  (5.26),  where  the  change  to  the  kth 


degree  of  freedom  for  the  ith  mode  in  the  absence  of  Coriolis  effects  is: 


P  f  n  ’(’kj  f 


If  Coriolis  effects  are  present,  then  Equation  (5.49)  can  be  used  to  compute  the  perturbed 


complex  eigenvectors. 


The  above  equation  is  a  linear  approximation  of  the  perturbations  in  the 


eigenvectors.  Using  the  results  of  reanalvsis,  one  obtains  the  full,  nonlinear  changes  in  the 


eigenvector.  These  two  candidate  procedures  have  definite  trade  offs.  Using  Equation 


(6.6)  avoids  another  finite  element  run.  However,  Equation  (6.6)  involves  many  matrix 


multiplications  and  manipulations  and  additional  programing.  Doing  the  rear,  aly sis 


involves  another  finite  element  run,  but  it  is  simple  to  do  and  provides  the  exact 


njK  n*  rkw.  www  tw a*  tv*  -ch  tv*  tjovti-jc  rjt  t vmrr.  jtvtk.  vt  vr*  yty^urr^  v»  ut  mt  ltvwti  lti  vt  vi 
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intermediate  answers.  It  also  incorporates  the  changes  in  centrifugal  force  due  to  mass 
changes,  which  is  something  Equation  (6.6)  cannot  do.  As  a  side  benefit,  the  intermediate 
result  of  the  predictor  method  alone  are  provided. 

A  third  approach  that  one  can  use  in  the  corrector  step  is  to  simply  use  the 
original,  unperturbed  eigenvector.  This  should  provide  reasonable  results  provided  that  the 
change  between  the  original  and  objective  systems  is  small.  The  benefit  of  this  procedure 
is  that  neither  a  finite  element  reanalysis  nor  a  new,  intricate  matrix  solution  is  required. 

As  shown  in  Figure  11,  the  finite  element  solver  is  first  run.  This  is  done  to 
generate  the  system  structural  matrices  and  the  eigenvectors.  Then,  using  this 
information,  the  equation  of  constraint  for  the  predictor  is  used  using  the  frequency  change 
equation  or  equations  given  in  Equation  (6.4)  or  (6.5)  and  (6.6).  This  constraint  is  used  in 
the  next  step,  the  optimizer.  Optimization  is  accomplished  with  respect  to  minimum 
weight  or  minimum  change.  Finally,  a  finite  element  analysis  is  done  to  obtain  the  system 
matrices  and  eigenvectors  for  the  perturbed  system  produced  by  the  predictor. 

Figure  12  shows  an  overview  of  the  corrector  step.  First,  the  finite  element  solver 
is  run  for  the  perturbed  system  resulting  from  the  predictor.  This  is  done  to  generate  the 
system  structural  matrices  and  the  eigenvectors.  Then,  using  this  information,  the 
equation  of  constraint  for  the  predictor  is  used  using  the  energy  balance  equation  or 
equations  given  in  Equation  (6.4)  or  (6.5)  and  (6.6).  This  constraint  is  used  in  the  next 
step,  the  optimizer.  Finally,  a  finite  element  analysis  is  done  to  obtain  the  system 
matrices  and  eigenvectors  for  the  perturbed  system  produced  by  the  predictor.  Note  that 
this  last  finite  element  analysis  is  not  really  required;  it  is  done  to  confirm  the  results. 

Only  the  initial  analysis  and  the  intermediate  analysis  are  required. 


Examples 

The  first  physical  problem  to  which  the  predictor-corrector  method  is  applied  is 
shown  in  Figure  10.  It  is  an  initially  uniform  aluminum  2024T-6  beam  of  length  L  equal 
to  250  mm,  with  moment  of  inertia  I  equal  to  3.2552E4  mm4,  cross  sectional  area  A  of 


i 


$ 
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Figure  12.  Corrector  Overview 


625  mm2,  Young’s  modulus  E  of  73.77E3  MPa,  Poisson’s  ratio  v  equal  to  0.33,  and 
density  p  of  2.774E-9  Mg/mm3.  The  beam  rotates  at  a  speed  fl  of  80  h z. 

A  preliminary  step  done  here  is  to  determine  the  fundamental  flexural  frequency 
for  the  problem  of  the  nonrotating  beam.  This  will  be  used  to  test  the  convergence  of  the 
finite  element  mesh.  The  exact  solution  for  fundamental  frequency  is  3.5043E3  rad/sec 
(see  Appendix  B).  Table  1  shows  that  the  accuracy  of  the  mesh  improves  when  more 
elements  are  used. 


Table  1.  Mesh  Convergence 


For  the  rotating  beam,  the  bending  effects  of  the  beam  are  separated  from  the 
axial  effects.  This  is  done  because  the  effect  of  extensional  (membrane)  stiffness  on  the 
predictor-corrector  equations  is  linear  in  the  design  variable  ae  while  the  effect  of  thickness 
on  bending  stiffness  is  cubic  in  ae.  The  separation  has  been  accomplished  in  an  artificial 
way  by  splitting  each  element  into  two  superimposed  subelements:  one  with  moment  of 
inertia  and  no  cross  sectional  area  and  another  with  cross  sectional  area  and  no  moment  of 
inertia.  At  first  this  sounds  unmotivated,  but  it  is  consistent  with  the  additive  value  of 
stiffness.  Bending  moment  of  inertia  and  cross-sectional  area  are  indeed  distinct 
properties.  Any  relationship  between  the  two  is  a  consequence  of  a  specific  geometry 
where  the  two  quantities  are  related  by  a  common  parameter.  For  a  rectangular  beam, 
both  the  moment  of  inertia  and  cross  sectional  area  are  dependent  on  thickness. 

Therefore,  in  this  case,  this  parameter  is  thickness,  and  optimization  of  this  parameter  will 
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be  done. 

The  beam  is  divided  into  four  equal-length  finite  elements,  and  each  element  is 
composed  of  two  subelements.  The  mesh  is  shown  in  Figure  13. 

For  the  given  mesh,  the  fundamental  frequency  of  the  rotating  beam  is  oj1  = 

3.3741 17E3  rad/sed.  The  modal  frequency,  equals  l.lwj  and  will  therefore  be 
3.711529E3  rad/sec  with  a  desired  change  of  frequency  of  3.374120E2  rad/sec.  The 
predictor  obtains  a  fundamental  modal  frequency  w \  of  3.722036E3  rad/sec.  This  is  within 
3.11%  of  the  desired  frequency  change.  The  corrector  step  is  then  done,  using  the  method 
that  employs  an  updated  eigenvector  from  the  intermediate  MSC/NASTRAN  analysis. 

The  fundamental  frequency  from  the  corrector,  uj,  is  3.713899E3  rad/sed.  Thus  it  is  seen 
that  the  predictor-corrector  obtains  the  frequency  change  to  within  1%  of  the  desired 
change. 


4 
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Figure  13.  Finite  Element  Model  for  Rotating  Beam 
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The  next  test  case  for  optimal  redesign  will  be  the  simple  rotating  blade  from 
Chapter  4.  The  blade  is  made  of  Inconel  718  steel  with  Young’s  modulus  of  202.696E3 
MPa,  Poisson’s  ratio  of  0.29,  and  density  of  8.2121E-9  Mg/mm3.  The  blade  outer  radius  is 
787.4  mm  as  measured  from  a  the  center  of  a  rigid  hub  of  radius  393.70  mm.  has  a  width 
of  556.78  mm,  and  is  initially  of  uniform  thickness  87.38  mm.  The  rotation  speed  0  is 
81.2  hz. 

Figure  14  shows  a  finite  element  model  for  this  blade,  which  is  the  same  as  the 
model  in  Figure  2.  The  rigid  hub  is  not  modeled.  This  finite  element  model  consists  of  four 
quadrilateral  elements,  and  each  element  is  once  again  broken  into  two  subelements.  The 
first  subelement  has  only  membrane  properties  and  the  second  one  has  only  bending 
properties.  As  in  the  case  of  the  beam,  this  is  done  in  order  to  separate  membrane 
stiffness  from  bending  stiffness. 

A  difficulty  presented  by  quadrilateral  elements  in  many  general  purpose  finite 
element  programs  is  the  that  stiffness  and  “coupled”  mass  matrices  output  are  in  element 
coordinates.  (The  “coupled”  mass  matrix  is  a  consistent  mass  matrix  with  the  rotational 
terms  removed.)  Since  the  eigenvectors  are  always  in  the  basic  (global)  coordinate  system, 
it  is  necessary  to  transfer  the  stiffness  and  mass  matrices  to  the  basic  system.  The 
stiffness  matrix  can  be  transformed  by: 

nw  *  mT[Ke]ement]m  (6.io) 

where  lKbasic]  is  the  element  stiffness  matrix  in  the  basic  system,  [Kelement]  is  the 
element  stiffness  matrix  in  the  element  system,  and  [T]  is  the  element-to-basic 
transformation  matrix.  The  element  mass  matrices  are  similarly  transformed  to  the  basic 
coordinate  system. 

The  fundamental  modal  frequencj'  for  the  rotating  blade  of  Figure  6.4.  u>1.  is 
2.478E3  rad/sec.  For  a  10%  desired  increase  in  ui1,  uf  is  2.726E3  rad/sec  and  A^  is 
2.483E2  rad/sec.  The  predictor  gives  of  2.653E3  rad/sec.  Au^  is  1.749E2  rad/sec, 
which  is  70.42%  of  the  desired  frequency  change.  The  corrector  obtains  u/j  of  2.705E3 
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Figure  14.  Finite  Element  Model  for  Rotating  Blade 


rad/sec.  is  2.268E2  rad/sec  which  is  91.33%  of  the  desired  change. 

Figure  15  shows  the  cross  section  of  the  optimal  stuctural  design  of  the  blade.  The 
original  structure  is  shown  in  solid  lines  and  the  optimal  structure  is  in  dashed 
lines.  Thickness  differences  have  been  exaggerated.  Refined  meshes  will  be  the  subject  of 
the  next  section. 

To  summarize  the  results  of  the  previous  two  examples,  the  governing  equation  of 
the  predictor  step,  Equation  (6.2),  is  a  linear  representation  of  a  nonlinear  perturbation  of 
the  mass  and  stiffness  matrices  and  the  matrix  of  eigenvectors  with  respect  to  a  change  in 
thickness.  However,  the  higher  order  perturbation  terms  are  important.  The  predictor 
achieved  only  70%  to  96%  of  the  desired  frequency  change.  Additional  nonlinearities  are 
brought  back  into  the  analysis  in  the  corrector  in  the  form  of  recalculated  eigenvectors. 
Changes  in  centrifugal  force  due  to  mass  changes  are  also  incorporated  in  the  reanalysis. 


Figure  15.  Optimized  Shape  of  Rotating  Blade 


The  predictor-corrector  method  outlined  in  this  section  obtained  90%  to  99%  of  the  desired 
frequency  change  in  one  complete  cycle.  The  end  result  is  an  improved  method  for  the 
optimal  redesign  of  a  rotating  structure. 


CHAPTER  7 


OPTIMAL  REDESIGN  OF  ADVANCED  SYSTEMS 

Review  of  Solution  Procedure 

A  method  has  been  devised  to  obtain  the  optimal  stuctural  redesign  of  a  rotating 
system  that  incorporates  centrifugal  and  Coriolis  effects.  This  method  extends  the  inverse 
perturbation  scheme  devised  by  Hoff  into  the  nonlinear  phj'sical  problem  regime.  We  must 
deal  with  the  structurally  nonlinear  effects  of  rotation  and  to  obtain  the  optimal  redesign 
with  a  minimum  number  of  finite  element  analyses.  In  order  to  account  for  the  effect  of 
the  static  displacements  due  to  rotation  on  the  bending  frequencies  and  mode  shapes,  a 
new  finite  element  solution  sequence  was  created.  This  involved  a  preliminary  static 
solution  where  the  displacements  due  to  rotation  were  first  calculated.  Then,  a  differential 
stiffness  matrix,  [KD],  was  calculated.  This  matrix  represents  the  effect  of  the 
displacements  on  the  structural  stiffness.  A  Coriolis  matrix  was  also  generated  that 
incorporated  the  energy-conserving  coupling  effects  of  velocity  in  one  plane  to  displacement 
in  another  plane.  The  differential  stiffness  matrix  and  the  Coriolis  matrix  were  then 
incorporated  into  the  modal  analysis  and  the  frequencies  and  mode  shapes  were  obtained. 

For  the  rotating  problem  without  Coriolis  effects,  Equat  <  n  (3.2)  gives  the  scalar 
relationship  between  the  element  change  ae  and  the  desired  change  in  natural  frequency 
for  the  ith  mode,  Awj.  This  equation  was  used  as  an  equality  constraint  for  the  function  to 
be  opimized,  Equation  (6.3).  The  perturbed  eigenvectors  for  the  new  system  were  then 
obtained  via  an  intermediate  finite  element  reanalysis.  These  were  used  in  Equation  (6.4) 
to  correct  the  energy  imbalance  between  the  predicted  system  and  the  desired  system 
through  additional  structural  changes.  This  equation  was  used  as  an  equality  constraint  in 


the  optimization  of  Equation  (6.3).  Finally,  the  new  eigenvalues  and  eigenvectors  of  the 
redesign  system  were  obtained  by  finite  element  reanalysis. 

Rotating  Compressor  Blade 

The  next  step  is  to  use  the  above  method  on  more  complicated  problems.  Figure  16 
shows  a  finite  element  model  for  a  curved  rotating  blade.  The  blade  is  made  of  Inconel  718 
steel,  has  a  radius  of  254.0  mm,  a  length  of  69.34  mm,  and  rotates  at  a  speed  of  200  hz. 

It  has  an  angle  of  attack  of  30  degrees  and  is  modeled  after  a  NACA  64  airfoil.  This  is  a 
blade  typically  found  in  a  jet  engine  high-pressure  compressor. 


z 


Figure  16.  Rotating  Blade 

The  finite  elements  are  each  divided  into  two  subelements;  one  with  membrane 
properties  only  and  another  superimposed  element  with  only  bending  properties.  This 
finite  element  model  has  approximately  1000  degrees  of  freedom.  The  elements  are 
grouped  (linked)  into  twelve  regions.  During  the  analysis,  the  thickness  of  the  regions  will 


***** 


k 


i  4'M’i  |'|j 


change,  but  the  elements  within  each  region  will  maintain  a  common  thickness.  (Each 
superimposed  membrane  and  bending  element  region  will  also  keep  a  common  thickness. 
The  twelve  regions  are  shown  in  Figure  17.)  Table  2  shows  the  initial  thickness,  in 
millimeters,  of  the  regions. 


Figure  17.  Blade  Regions 

Let  Qj  be  the  element  change  property  for  the  ith  element.  Using  Equation  (6.2)  for 
the  first  mode,  define  the  contribution  of  the  ith  element  Ej  to  the  frequency  change 


E1  =  M7  +  "W  * 


-  w/WjEmj]  {tl>h 


Now  let  be  the  common  element  change  property  for  all  of  the  elements  in  region 


> 
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Table  2.  Initial  Region  Thickness 


Region  Number 

Thickness  (mm) 

1 

1.734 

2 

1.734 

3 

1.734 

4 

2.312 

5 

2.312 

6 

2.312 

7 

3.005 

8 

3.005 

9 

3.005 

10 

3.467 

11 

3.467 

12 

3.467 

k.  In  the  corrector,  the  lower  limit  on  ak  must  be  changed  so  that  no  region  will  have  a 
thickness  less  than  one-half  of  the  thickness  of  the  region  in  the  original  geometry.  Let  ak| 
denote  this  new  lower  bound  of  the  element  change  parameter  for  each  -egion.  Let  Ek  be 
the  sum  of  the  Ek  factors  for  each  element  in  region  k.  Therefore,  using  Equation  (6.2). 
the  optimization  problem  for  the  predictor  becomes  the  minimization  of: 

12 

f(ak>=X>k)2  (7-2) 

k=  1 

subject  to  the  equality  constraint: 

12 

£Ekak  =  M  (7.3) 

k  =  1 


and  the  inequality  constraints: 
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ak  >  -0.5  (7.4) 

4  <  1E5  (7.5) 

For  the  corrector,  Equation  (4.4)  is  modified  as  follows.  Let  Jj  be: 

Ji  =  Mltk'iHn  (7.6) 

Define  G  to  be: 

G  =  w’fty'tftmHV-'}!  -  {V>,)T[k]{^,}1  (7.7) 


where  {i/i'}  is  the  perturbed  eigenvector,  «/'  is  the  desired  eigenvalue,  and  [k]  and  [m]  are 
the  global  stiffness  and  mass  matrices,  respectively. 

Now  let  Qk  be  the  common  element  change  property  for  all  of  the  elements  in  region 
k.  Let  Jk  be  the  sum  of  the  Jk  factors  for  each  element  in  region  k.  Therefore,  using 
Equation  (6.5),  the  optimization  problem  for  the  predictor  becomes  the  minimization  of: 


12 

f(ak)  =  £(ak)2  (7.8) 

k=l 

subject  to  the  equality  constraint: 

12 

EJk°k  =  G  (7.9) 

k=  1 

and  the  inequality  constraints: 

“k  >  Qk[  (7-10) 

Qk  <  1E5  (7.11) 


For  this  problem  another  FORTRAN  postprocessor  must  be  written  that  groups 
together  the  stiffness,  mass,  differential  stiffness,  eigenvector,  and  transformation 
matrices  for  each  element.  This  is  done  so  as  to  minimize  matrix  storage  requirements 
and  programming  steps  during  the  predictor  or  corrector.  Rather  than  read  in  all  of  the 
matrices  for  all  of  the  elements  at  once,  they  are  read  in  per  element.  The  matrix 
operations  for  the  predictor  or  corrector  are  then  performed,  the  results  are  stored,  and 
then  the  program  loops  back  to  read  in  the  matrices  for  the  next  element. 

Three  complete  finite  element  static  and  dynamic  analyses  are  required  for  each 
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cycle  in  this  method.  The  first  run  gives  the  eigenvalues  and  mode  shapes  for  the  original 


system.  The  second  one  gives  the  values  for  the  system  output  by  the  predictor.  The  final 


analysis  verifies  the  results  of  the  corrector. 


Another  possible  objective  function  is  minimum  weight.  In  this  case,  the  function  to 


be  minimized  in  the  predictor  comes: 


fXtelJ)  -  2^Akak 

k=  1 


(7.12) 


subject  to  the  inequality  constraints: 


ak  >  -0.5 


(7.13) 


<  1E5 


(7.14) 


In  the  corrector  step  for  the  minimum  weight  problem,  the  objective  function 


remains  minimum  minimum  weight.  However,  as  in  the  minimum  change  case,  the  lower 


limit  on  the  element  change  parameter  must  be  changed  so  that  no  region  will  have  a 


thickness  less  than  one-half  the  thickness  of  the  region  in  the  original  geometry. 


Results  of  Analysis 


It  is  first  noted  that  for  the  nonrotating  system,  the  fundamental  frequency  is 


7.6656E3  rad/sec.  For  the  rotating  system  including  centrifugal  effects,  the  fundamental 


frequency  is  8.3802E3  rad/sec.  This  implies  a  9.32%  increase  in  fundamental  frequency 


due  to  the  centrifugal  effect  of  rotation. 


Figures  IS  and  19  show  the  first  two  eigenvectors  for  the  rotating  blade  when 


centrifugal  effects  are  not  included.  The  first  eigenvector  is  a  flapping,  or  bending,  mode 


and  the  second  one  is  a  twisting,  or  torsional  mode.  The  fundamental  eigenvector  is  indeed 


that  of  out-of-plane  bending.  Figures  20  and  21  show  the  first  two  eigenvectors  when 


centrifugal  effects  are  included.  Notice  that  the  mode  shapes  remain  unchanged. 


Therefore,  the  effect  of  rotation  changes  the  eigenvalues  but  the  eigenvectors  retain  their 


identity. 


Figure  18.  Mode  Shape  1,  Blade 
No  Rotational  Effects 


The  problem  of  the  rotating  blade  was  analyzed  for  several  cases  and  using  several 
different  methods  to  account  for  nonlinearities.  In  cases  1  and  2,  a  10%  increase  in  the 
fundamental  eigenfrequency  was  desired,  with  the  objective  function  for  case  1  being 
minimum  change  and  the  objective  function  for  case  2  being  minimum  weight.  In  cases  1 
and  2,  centrifugal  effects  were  included  in  both  the  structural  analysis  and  the 
optimization,  but  Coriolis  effects  were  neglected.  Results  of  both  predictor  and  corrector 
are  shown  in  Table  3.  Note  that  the  predictor  results  can  be  considered  to  be  results  from 
a  linear,  one-step  analysis  since  the  effect  of  redesign  on  the  eigenvectors  does  not  enter 
into  the  predictor  procedure.  Improvements  from  the  predictor  to  corrector  step  show  the 
benefit  of  the  use  of  the  nonlinear  optimization  techniques. 

In  Table  3,  the  first  colum  denotes  the  objective  functions  used  in  the  predictor  and 


Figure  19.  Mcde  Shape  2,  Blade 
No  Rotational  Effects 


the  corrector  respectively.  The  symbol  C/C  denotes  minimum  change  in  both  parts.  If 
W/W  is  indicated,  minimum  weight  was  used  in  both  the  predictor  and  corrector.  The  use 
of  W/C  which  symbolizes  that  minimum  weight  was  used  in  the  predictor  while  minimum 
change  was  used  in  the  corrector  indicates  a  hybrid  system  (Case  3),  which  shall  be 
described  later.  The  next  column  entry,  Wj,  indicates  the  fundamental  frequency  for  the 
system.  The  desired  frequency  is  denoted  by  oij.  The  fundamental  frequency  resulting 
from  the  predictor  geometry  is  indicated  by  The  percentage  of  the  desired  frequency 
change  accomplished  by  the  predictor  is  symbolized  by  A^.  The  percent  weight  change 
resulting  from  the  predictor,  AWp,  is  given  in  the  next  column.  The  fundamental 
frequency  of  the  corrector  geometry  is  denoted  by  J\.  The  percentage  of  the  desired 
frequency  change  accomplished  by  the  corrector  is  given  by  %  Ac.  In  the  final  column,  the 


Figure  20.  Mode  Shape  1,  Blade 
Rotational  Effects  Included 

percent  weight  change  resulting  from  the  corrector,  %AWC,  is  given. 

Figure  22  shows  the  final  spanwise  optimized  thickness  of  the  structure  for  Case  1. 
Figure  23  shows  the  final  optimized  shape  of  the  structure  for  Case  2.  Notice  that  in  the 
minimum  change  example,  emphasis  is  given  to  adding  material  at  the  root.  In  the 
minimum  weight-minimum  weight  example  (Case  2),  all  of  the  regions  except  for  the  root 
have  been  reduced  to  the  lower  limit  on  thickness.  This  is  the  pathological  case  in 
optimization  where  the  system  is  driven  to  an  extreme.  When  this  is  done  in  this  example, 
undesirable  side  effects  occur,  such  as  mode  switching.  The  first  bending  mode  is  no  longer 
the  fundamental  frequency  and  the  the  solution  to  the  problem  in  optimization  is  no  longer 
dependable.  The  frequency  results  shown  for  the  corrector  are  for  the  bending  mode: 
however,  this  frequency  is  technically  no  longer  u>j.  A  way  out  of  this  quandary  can  be 
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Figure  21.  Mode  Shape  2,  Blade 
Rotational  Effects  Included 


found  by  a  close  examination  of  the  results  of  the  predictor.  This  step  obtains  99%  of  the 
desired  frequency  change  and  also  a  weight  reduction  of  —14%.  Therefore,  this  solution  is 
close  to  the  frequency  constraint,  and  only  a  small  change  is  necessary  to  satisfy  it.  This 
implies  that  a  hybrid  approach  involving  a  predictor  stop  with  a  minimum  weight  objective 
function  and  a  corrector  with  a  minimum  change  objective  function  could  work. 

The  results  of  this  hybrid  step  are  shown  in  Case  3  and  Figure  24.  In  this  example, 
material  is  added  at  the  root  but  proportionally  less  than  in  the  minimum  change- 
minimum  change  situation.  Emphasis  is  given  to  removing  material  from  the  outboard 
regions,  with  most  material  removed  from  the  second  set  of  elements  from  the  end.  Since 
minimum  weight  is  the  objective  of  the  predictor,  it  is  not  suprising  that  more  material  is 
removed  in  Case  3  than  in  the  Case  1  situation. 


Cases  involving  degradations  of  the  solution  procedure  are  shown  in  Appendix  C. 
These  were  done  to  observe  the  effect  of  the  absence  of  centrifugal  effects  on  the  problem 
in  optimization. 

Two  other  problems  were  studied;  both  involved  large  (30%)  changes  in  the 
fundamental  frequency.  In  Case  4,  the  30%  change  is  accomplished  in  one  step.  A  second 
iteration  is  performed  to  obtain  an  improved  solution  (Case  5).  In  another  situation,  the 
30%  change  is  broken  down  into  three  10%  increments  (Cases  6  through  8).  In  both  of 
these  examples,  only  a  minimum  change  optimization  function  is  used.  Table  4  shows  the 
results  of  the  iterative  procedure.  The  linear  predictor  step  obtains  the  desired  frequency 
change  with  less  than  24%  accuracy,  but  at  the  end  of  the  first  iteration,  the  desired 
frequency  change  is  accomplished  to  within  less  than  1%.  The  second  iteration  is  done  for 
completeness,  and  gives  the  desired  change  in  fundamental  frequency  to  within  — of  1%. 

In  Table  5,  each  iteration  obtains  the  desired  change  in  frequency  for  that  iteration 
to  less  than  1%.  The  final  iteration,  which  completes  the  30%  change,  gives  the  desired 
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change  to  within  -j-QQ  of  1%.  These  two  Tables  show  that  excellent  accuracy  on  the 
frequency  goal  is  obtained,  showing  the  feasibility  of  making  large  changes. 

The  use  of  the  word  “accuracy”  needs  some  close  examination.  “Accuracy”  has 
been  used  to  describe  how  close  the  frequency  change  obtained  by  solving  the  problem  in 
optimization  is  to  the  desired  frequency  change.  This  is  done  in  terms  of  a  discrete  finite 
element  model.  An  analysis  was  made  in  Appendix  B  of  the  convergence  of  discrete  finite 
element  models  to  a  theoretical  solution  involving  a  nonrotating  beam.  This  can  be  used  to 
give  an  indication  of  the  accuracy  of  the  finite  element  model  in  that  case.  However,  no 
such  comparison  to  a  theoretical  solution  is  made  for  the  rotating  blade  presented  in  this 
chapter.  Closed-form  solutions  to  non-abstract  structural  systems  are  difficult  to  obtain. 
Therefore,  in  all  discussions  on  accuracy  in  this  chapter,  comparisons  are  made  from  one 


discrete  model  to  another. 
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Figure  22.  Optimally  Redesigned  Blade,  Case  1 
Minimum  Change 


Optimization  Constraints 


The  optimizer  is  run  iteratively,  with  the  solution  of  the  previous  step  becoming  the 


starting  point  for  the  succeeding  step  (the  initial  starting  point  is  the  origin).  In  nonlinear 


mathematical  programming,  the  approach  is  to  minimize  the  objective  function  while 


driving  the  equality  constraint  function  to  zero.  Table  6  shows  for  each  case  the  number  of 


iterations  i  required  for  final  solution,  the  value  of  the  objective  function  f,  and  the  value  of 


the  constraint  A.  The  solution  is  terminated  when  further  iterations  obtain  improvements 


only  of  less  than  0.001  (Reason  1)  or  further  iterations  cause  the  solution  to  diverge,  as 


indicated  by  increasing  values  of  objective  function  or  constraint  (Reason  2).  The  table  also 


gives  the  reason  for  termination,  1  or  2,  as  indicated  above.  Note  that  in  no  case  were 


more  than  seven  iterations  required. 


Figure  23.  Optimally  Redesigned  Blade,  Case  2 
Minimum  Weight 


In  Table  6,  it  is  noticed  that  cases  with  a  minimum  change  objective  function  have 
the  optimized  value  of  the  function  A  close  to  zero,  but  cases  with  a  minimum  weight 
objective  function  often  have  negative  values.  This  is  logical:  for  minimum  change,  low 
values  imply  low  sum  of  the  squares  of  change.  For  minimum  weight,  large  negative 
values  imply  a  large  reduction  of  volume. 


Rotating  Beam  Incorporating  Coriolis  Effects 
The  next  example  will  consider  the  case  of  the  rotating  beam  presented  in  the 
prevous  chapter  (Figure  13)  at  high  rotational  speed.  Once  again,  the  beam  will  be  divided 
into  four  elements  with  each  element  composed  of  two  subelements:  one  having  bending 


properties  only  and  the  other  posessing  only  tension-compression  properties.  The  rotation 
speed  will  be  300  hz. 


Figure  24.  Optimally  Redesigned  Blade,  Case  3 
Hybrid 

For  this  problem,  the  fundamental  eigenfrequency  without  rotation  is  2.093729E3 
rad/sec.  Figure  25  illustrates  the  first  mode  shape  of  the  beam,  a  bending  mode.  When 
centrifugal  effects  are  included,  the  fundamental  frequncy  is  2.952001E3  rad/sec,  an 
increase  of  40.99%.  When  Coriolis  effects  are  included  in  addition,  the  fundamental 
frequency  drops  slightly  to  2.951942E3  rad/sec.  The  inclusion  of  Coriolis  forces  in  the 
rotating  problem  decreases  the  fundamental  frequency  by  -2.00E-3%  from  the  rotating 
problem  that  includes  centrifugal  but  not  Coriolis  effects. 

Table  8  summarizes  the  results  from  the  optimal  redesign  of  the  rotating  beam. 
The  notation  is  the  same  as  for  the  rotating  blade.  In  Case  9,  centrifugal  effects  are 
included  in  both  the  structural  analysis  and  in  the  optimization.  A  minimum  change 
objective  was  used.  Case  10  was  identical  to  case  1;  however,  the  hybrid  case  was 


Figure  25.  Mode  Shape  1,  Beam 

analyzed  incorporating  a  minimum  weight  objective  function  in  the  predictor  and  a 
minimum  change  objective  function  in  the  corrector.  Coriolis  effects  were  included  in  Case 
1 1  in  both  the  structural  analysis  and  equations  of  constraint.  Optimization  was 
accomplished  using  a  minimum  change  objective  function.  Figures  26  and  27  show  the 
final  configuration  of  the  rotating  beam  for  Cases  9  aiid  10,  respectively.  Case  11.  being  a 
minimum  change  case,  shows  an  optimized  shape  similar  to  Figure  26. 

Optimization 

As  in  the  previous  section,  the  values  for  constraint  functions,  constraint  values,  and 
reason  for  termination  are  given  for  each  of  the  three  cases  in  Table  8.  For  the  cases 
involving  complex  equations,  the  constraint  shown  is  the  real  constraint. 

Notice  that  the  function  evaluations  for  the  minimum  change  objective  are  close  to 
zero,  and  the  function  evaluations  for  the  minimum  weight  objective  (predictor  only)  are 
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Figure  26.  Optimally  Redesigned  Beam,  Case  9 
negative.  They  are  not  as  big  in  this  example  as  in  the  blade  example. 

Summary  of  Results 

The  predictor-corrector  methods  breaks  the  solution  of  the  problem  of  nonlinear 
optimal  redesign  into  two  parts.  The  first  part,  the  predictor,  solves  for  the  required 
structural  changes  for  a  given  required  change  of  frequency.  In  this  step,  the  effect  of 
structural  changes  on  the  mode  shapes  is  not  considered.  Therefore,  this  part  of  the 
solution  may  be  considered  as  a  conventional  linear  structural  analysis.  In  the  corrector, 
the  effect  of  the  structural  changes  on  the  mode  shapes  is  taken  into  account  and  the 
system  is  once  again  modified  to  obtain  an  improved  solution. 

In  all  of  the  examples  involving  the  rotating  blade,  the  final  result  of  the  predictor- 
corrector  approacn  obtains  the  desired  frequency  change  to  within  one  percent.  In  the 
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Figure  27.  Optimally  Redesigned  Beam,  Case  10 


minimum  change  cases,  the  linear  predictor  overshoots  the  solution  by  a  few  percent.  The 
corrector  changes  the  final  solution  so  that  the  eigenfrequency  is  at  the  desired  value. 

Even  for  large  changes,  the  predictor-corrector  method  obtains  the  desired  solution,  if 
suitable  iteration  or  incrementing  is  done. 

For  the  rotating  beam  with  a  minimum  weight  objective  function,  the  linear  predictor 
undershoots  the  desired  frequency  goal  by  quite  a  bit.  as  much  as  35*70.  The  corrector 
obtains  the  desired  frequency  to  within  one  percent. 

The  rotating  beam  shows  some  other  interesting  results.  In  Case  9,  the  desired 
change  is  obtained  to  within  4%.  In  Case  10,  there  is  a  lot  of  undershoot  by  the  predictor, 
but  the  corrector  obtains  the  desired  solution  to  within  6%.  However,  when  both  centrifugal 
and  Coriolis  forces  are  included  in  Case  11,  the  best  solution  is  found.  The  linear  approach 
gives  an  answer  to  within  5%  and  the  corrector  improves  this  to  within  1%.  The  method 
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free  end  has  the  minimum  thickness  and  the  end  bulges  out,  though  it  remains  less  thick 
than  the  original  design.  This  was  also  obtained  by  Queue  and  Trompette  (Figure  28), 


In  work  on  the  nonrotating  beam,  Karihaloo  and  Niordson  (1973)  obtained  the  classical 
tapered  beam  for  the  optimum  shape  of  the  initially  uniform  beam,  where  the  first  mode  is 


of  interest.  This  was  confirmed  by  OlhofT  and  Parbery  (1976)  and  expanded  to  higher 
order  modes  (Olhcff,  1977). 


used  in  Case  11  represents  the  best  theoretical  formulation.  The  equations  used  represent 
the  full  nonlinear  structural  approach  with  both  centrifugal  and  Coriolis  effects. 

One  of  the  reasons  behind  the  lack  of  accuracy  surrounding  the  problem  of 
optimization  with  respect  to  minimum  weight  is  that  this  optimization  function  in  the 
presence  of  the  nonlinear  effects  of  the  corrector  drives  most  of  the  structure  to  the 
pathological  extreme  of  the  lower  bound  on  thickness.  This  is  an  undesirable  solution.  Not 
only  does  it  introduce  inaccuracies  that  somewhat  distance  the  frequency  change  from  the 


desired  value,  but  it  also  causes  mode  switching  that  results  in  the  bending  mode  no  longer 
being  the  fundamental  mode.  Since  higher  order  modes  were  not  tracked  in  this  procedure, 
the  solution  to  the  problem  in  optimization  is  no  longer  valid. 

Comparison  with  Other  Methods 

Queau  and  Trompette  (1981)  obtained  minimum  weight  designs  with  constraints  on 
frequency.  Their  method  incorporated  the  centrifugal  effects  but  not  Coriolis.  In  the 
method  implemented  here,  when  minimum  weight  is  employed,  the  second  station  from  the 


with  the  dashed  line  indicating  the  final  optimized  shape. 

Olhoff  and  Parbery  (1982)  examined  the  optimization  of  rotating  beams  with  respect 
to  frequency  constraints.  However,  they  employed  lumped  masses  which  tend  to  alter  the 
optimized  shape  from  the  pureljr  distributed  mass  approach.  Their  final  shapes  indicated 
tapering  except  near  the  region  surrounding  a  lumped  mass  where  bulging  then  occured. 
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CHAPTER  8 


CONCLUSIONS 

General  Conclusions 

The  effect  of  rotation  on  the  eigenvalues  and  eigenvectors  of  a  system  was 
successfully  modeled  through  the  use  of  finite  element  theory.  Both  centrifugal  forces  and 
Coriolis  effects  were  included.  The  centrifugal  forces  in  this  work  were  treated  as  a  static 
load  and  therefore  the  effect  was  included  as  a  mild  geometric  nonlinearity  in  the 
differential  stiffness  matrix.  The  velocity-dependent  Coriolis  terms  were  formulated  into  a 
separate  Coriolis  matrix.  These  terms  were  found  to  be  very  small  in  the  examples 
treated  here.  The  results  obtained  by  the  finite  element  formulation  were  shown  to  be  in 
close  agreement  to  those  obtained  by  classical  means. 

By  far,  the  major  rotary  phenomenon  affecting  the  fundamental  eigenvalue  was 
centrifugal  force.  The  centrifugal  effect  increased  the  fundamental  frequency  of  the  blade 
and  the  beam  by  several  percent.  The  mode  shapes  remained  unchanged  from  the 
nonrotating  system.  The  centrifugal  forces  serve  to  stiffen  the  rotating  body.  Coriolis 
forces  slightly  decreased  the  fundamental  frequency.  Since  Coriolis  forces  are  velocity 
dependent,  they  affected  both  the  amplitude  and  phase  of  the  eigenvector  components 
creating  a  complex  eigenvector  problem.  However,  the  eigenvalues  remained  real. 

The  predictor-corrector  method  for  structural  optimization  using  inverse 
perturbation  was  extended  to  incorporate  centrifugal  forces.  To  accomplish  this,  the 
element  stiffness  matrices  were  separated  to  isolate  the  effects  of  in-plane  stiffness 
(membrane  or  extensions!),  bending  stiffness,  and  differential  stiffness.  When  the 
frequency  constraint  involved  a  ten  percent  increase  in  the  fundamental  frequency,  the 
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linear  predictor  obtained  the  required  change  within  two  percent  for  the  blade  or  five 
percent  for  the  beam.  The  nonlinear  corrector  obtained  a  final  optimized  system  that  met 
the  frequency  constraint  within  one  percent.  Thus  the  predictor-corrector  method  for 
nonlinear  redesign  obtained  excellent  agreement  between  the  desired  eigenvalue  and  the 
calculated  eigenvalue,  when  centrifugal  effects  were  included. 

When  Coriolis  effects  were  included,  both  the  magnitude  and  phase  of  the 
components  of  the  fundamental  eigenvector  were  required  to  obtain  the  equations  of 
constraint.  This  problem  in  complex  eigenvalue  analysis  was  also  adapted  to  the  nonlinear 
inverse  perturbation  predictor-corrector  approach.  The  method  was  applied  to  the  problem 
of  the  rotating  beam.  Once  again,  the  desired  frequency  change  was  obtained  to  within 
one  percent. 

The  problem  of  large  frequency  change  (30%)  was  tried  for  the  rotating  blade 
incorporating  centrifugal  effects.  Both  an  iterative  and  an  incremental  solution  were 
accomplished.  The  iterative  approach  obtained  the  desired  frequency  almost  exactly  with 
only  two  iterations  required.  The  incremental  approach  also  achieved  the  desired 
frequency  using  three  incremental  steps.  Thus  it  is  seen  that  the  predictor-corrector 
method  is  extraordinarily  stable,  obtaining  even  large  changes  with  excellent  correlation 
between  the  desired  change  in  the  eigenvalue  and  the  calculated  change  resulting  from  the 
redesign  process. 

Good  solutions  resulting  in  excellent  agreement  between  the  desired  frequency 
change  and  the  actual  frequency  change  always  obtain  when  a  minimum  change  objective 
function  is  used  in  both  the  predictor  and  the  corrector.  Minimum  weight  is  a  useful 
objective  function  for  the  predictor,  but  when  it  is  used  in  the  corrector,  a  pathological 
solution  results  with  all  mass  concentrated  at  one  area.  This  is  not  acceptable.  Not  only  is 
this  an  unrealistic  redesign,  but  the  bending  frequency  drifts  away  from  the  desired  value 
and  mode  switching  may  result.  To  correct  this  deficiency  in  the  minimum  weight 
approach,  a  minimum  change  objective  function  was  used  in  the  corrector  step.  This  hybrid 


approach  combines  the  desired  goal  of  minimum  weight  with  the  stability  of  the  minimum 
change  objective  function. 

The  solutions  calculated  by  the  predictor -corrector  method  were  compared  to  those 
obtained  by  other  authors.  These  solutions  were  for  minimum  weight  objective  under  a 
variety  of  frequency  goals.  Good  correlation  between  the  predictor-corrector  method  and 
other  methods  was  observed. 

In  summary,  the  predictor-corrector  method  for  optimal  redesign  as  extended  in  this 
dissertation  obtained  the  desired  frequency  changes  with  excellent  accuracy.  MSC/ 
NASTRAN  finite  element  solution  sequences  were  modified  to  incorporate  rotary  effects. 
DMAP  and  FORTRAN  languages  were  used  to  implement  the  derived  equations  of 
constraint.  Automated  Design  Synthesis  (ADS)  was  used  to  obtain  the  optimum  solution. 
The  methods  used  were  applied  to  several  test  problems,  one  being  a  curved  blade  with 
nearly  one  thousand  degrees  of  freedom.  In  each  case,  the  desired  frequency  change  was 
obtained  to  within  a  few  percent.  Therefore,  the  approach  works  and  can  be  applied  to  a 
variety  of  problems  in  optimal  redesign. 

Dissertation  Contribution 

The  major  contribution  is  the  extension  of  the  inverse  perturbation  method  to  include 
both  centrifugal  and  Coriolis  effects  in  the  optimization  of  rotating  systems  under 
frequency  constraints.  Either  minimum  structural  change  or  minimum  weight  objective 
functions  can  be  employed.  Excellent  agreement  between  the  desired  change  in  the 
fundamental  natural  frequency  and  the  actual  change  as  determined  by  reanalysis  was 
obtained  in  the  test  cases  tried.  Previous  researchers  have  neglected  the  Coriolis  effects  in 
the  optimization  of  rotating  systems.  Both  magnitude  and  phase-dependent  constraint 
equations  were  used  and  it  was  shown  how  the  equations  reduced  to  the  conventional 
rotating  solution  in  the  absence  of  Coriolis  effects. 

Another  original  contribution  is  the  implementation  of  a  method  for  modifying  an 
existing  large  scale  finite  element  program  to  incorporate  static  stiffening  loads  as 
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observed  in  the  rotating  frame  due  to  the  centrifugal  effect.  In  addition,  a  method  for 
incorporating  a  global  skew-symmetric  Coriolis  matrix  in  the  complex  modal  analysis  is 
developed.  Together,  these  provide  the  eigensolution  for  a  rotating  body.  The  technology 
developed  is  transportable  to  other  sites  by  anyone  familiar  with  MSC/NASTRAN  and 
DMAP. 

Suggestions  for  Future  Work 
Theoretical  Recommendations 

Coriolis  matrices  for  the  beam  element  using  a  consistent  mass  approach  and  in  a 
generalized  coordinate  system  should  be  obtained,  as  should  Coriolis  matrices  for  other 
elements.  The  problem  of  maximizing  the  differences  between  eigenvalues,  as  identified  by 
Olhoff  and  Parbery  (1982),  should  be  analyzed.  This  is  important  in  problems  where 
several  eigenvalues  are  tracked  in  the  redesign  process  with  emphasis  given  to  the 
difference  between  the  eigenvalues. 

Practical  Recommendations 

The  finite  element  method  employed  in  this  work  should  be  generalized  to  include 
more  types  of  elements,  such  as  higher-order  plates.  A  general  problem  incorporating 
several  different  elements  should  be  tried.  This  would  permit  the  analysis  of  a  real-world, 
built  up  structure  such  as  a  helicopter  blade.  An  improved  method  of  incorporating 
element  Coriolis  matrices  into  the  global  problem  within  the  context  of  existing  finite 
element  programs  should  be  devised.  This  should  be  done  so  that  the  employment  of  a 
Coriolis  matrix  in  complex  dynamic  analysis  would  be  transparent  to  the  user  and  would 
permit  the  automated  analysis  of  the  dynamic  problem  including  all  rotational  effects. 
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APPENDICES 


APPENDIX  A 


COMPUTER  IMPLEMENTATION 

Implementation  of  the  predictor-corrector  method  to  rotational  problems  is 
accomplished  by  a  combination  of  the  MSC/NASTRAN  program,  the  ADS  program  and 
FORTRAN  programming.  MSC/NASTRAN  is  used  to  accomplish  structural  analysis,  both 
static  and  modal,  and  generates  the  structural  matrices  (McNeal  1984).  FORTRAN 
programs  are  used  to  transform  the  structural  matrices  from  MSC/NASTRAN  into  a 
format  usable  for  general  analysis.  Then,  Direct  Matrix  Abstraction  Programming 
(DMAP)  is  used  to  generate  the  equations  of  equality  constraint.  Optimization  is  done  by 
use  of  the  ADS  program  (Vanderplaats,  1983). 

Figure  11  shows  the  flow  diagram  for  the  predictor  step.  The  finite  element  analysis 
is  accomplished  in  two  steps.  First  the  static  solution  of  MSC/NASTRAN.  SOL  24,  is  run 
to  generate  the  differential  stiffness  matrix  (MacNeal.  1981).  This  matrix  is  then 
checkpointed  to  the  restart  tape.  Then  SOL  3  is  run  for  the  real  modal  analysis  case  (or. 
alternatively,  SOL  28  is  run  for  the  complex  eigenvalue  problem).  The  problem  is  run  as  a 
restart  and  the  differential  stiffness  matrix  is  read  from  the  “old  problem  tape”  (Joseph, 
1984).  In  both  the  static  and  dynamic  runs,  the  mass,  stiffness,  differential  stiffness,  and 
eigenvector  matrices  are  written  on  tape  in  FORTRAN  compatible  binary  via  OUTPUT4. 
For  the  complex  run,  the  Coriolis  matrix  is  also  stored.  However,  MSC/NASTRAN  does 
not  itself  have  Coriolis  factors.  These  must  be  generated  by  a  FORTRAN  program.  The 
Coriolis  matrix  is  output  by  OUTPUT4  so  that  it  can  be  processed  by  the  same  programs 
as  are  the  other  matrices.  Finally,  a  nonrotating  dynamic  solution  is  done.  This  is  done  to 
generate  the  nonrotating  eigenvalues  and  also  to  produce  the  global  to  basic  system 
transformation  matrices  on  an  OUTPUT2  tape.  These  matrices  are  located  in  the  KDICT 
table. 

The  following  DMAP  ALTER  is  used  in  SOL  24  to  generate  the  differential  stiffness 
matrix.  Note  that  first  the  static  displacment  caused  by  the  centrifugal  forces  is  obtained 
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and  then  this  is  fed  back  into  the  EMG  processor  to  get  the  element  differential  stiffness 

matrices  (Gockel,  1984).  OUTPUT4  is  then  used  to  put  on  tape  the  matrix  of  the  element 

stiffness,  differential  stiffness,  and  mass  matrices. 

SOL  24  SSTATIC  ANALYSIS 
CHKPNT  YES 
DIAG  4,8,14 
ALTER  33 

OUTPUT 4  KELM.MELM//0/8  $ 

TAB PT  SIL/7  $ 

MATPRN  CSTM//  $ 

ALTER  30 
PRTPARM  III!  1  $ 

ALTER  97 
PRTPARM  III  1 1  $ 

ALTER  160 

EMG  EST, CSTM, MPT.DIT.,UGV.ETT.EDT/KDELM.KDDICT„.,/ 1/0/0/ 

///- 1/-  1/////////K6ROT  $ 

CHKPNT  KDELM.KDDICT  $ 

OUTPUT 4  KDELM//0/8  $ 

TABPT  KDDICT.,,,//  $ 

EMA  GPECT.KDDICT,KDELM,BGPDT,SIL.CSTM/KDNN./ -  1  $ 

CHKPNT  KDNN  $ 

EQUIV  KDNN.KDFF/SINGLE  $ 

CHKPNT  KDFF  $ 

COND  LBL3D,SINGLE  $ 

SCE1  USET.KDNN,„/KDFF,KDFS„„  $ 

CHKPNT  KDFF.KDFS  $ 

LABEL  LBL3D  $ 

EQUIV  KDFF,KDAA/OMIT  $ 

CHKPNT  KDAA  $ 

COND  LBL5D.OMIT  $ 

SMP2  USET.GO.KDFF/KDAA  $ 

CHKPNT  KDAA  $ 

LABEL  LBL5D  $ 

ADD  KDAA,KAA/KTOT/C,N.(1.0,0.0)/C,N,(1.0,0.0)  $ 

CHKPNT  KTOT  $ 

END  ALTER 
CEND 


This  following  DMAP  is  used  to  modify  SOL  3  to  produce  the  real  modal  solution  for 
the  rotating  problem  incorporating  differential  stiffness,  The  matrix  of  eigenvectors 
created  is  written  on  tape  by  OUTPUT4. 


SOL  3  $REAL  MODAL  ANALYSIS 
DIAG  4.8,14 
ALTER  395 
PRTPARM  lllll  $ 

ADD  KDAA,KAA/KNEW/C,N,(1.0,0.0)/C.N,(  1.0,0.0)  $ 
EQUIV  KNEW, KAA/ ALWAYS  $ 
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ALTER  448 

0UTPUT4  UGV,„,//0/8  $ 

END ALTER 
CEND 

For  the  complex  case  with  Coriolis  effects,  the  following  is  used  to  alter  SOL  28. 

The  global  Coriolis  matrix  is  read  via  INPUTT4.  Once  again,  OUTPUT4  is  used  to  write 

the  eigenvectors  on  tape. 

SOL  28  SCOMPLEX  MODAL  ANALYSIS 
DIAG  4,8,14 
ALTER  66 

INPUTT 4  /BCOR,,,,/ 1/8/0  $ 

MATPRN  BCOR/'/  $ 

EQUIV  BCOR.BGGX/ ALWAYS  $ 

SETVAL  //V,N,NOBGGX/0  $ 

ALTER  392 
PRTPARM  lllll  $ 

ADD  KDAA.KAA/KNEW/C,N,(1.0,0.0)/C,N,(1.0,0.0)  $ 

EQUIV  KNEW, KAA/ ALWAYS  $ 

ALTER  409 

MATPRN  KDXX,BDXX.MDXX//  $ 

ALTER  456 

OUTPUT4  UGV,,,,//-  1/8  $ 

MATPRN  UGV//  $ 

END  ALTER 
CEND 

In  applying  the  method  described  above  in  MSC/NASTRAN,  it  is  first  necessary  to 
divide  the  structural  elements  in  the  finite  element  model  into  two  elements:  one  with  only 
in-plane  (membrane  or  rodlike)  stiffness  and  one  with  only  out-of-plane  (bending)  stiffness. 
These  elements  are  then  superimposed.  The  element  stiffness  matrices,  created  by  DMAP 
module  EMG,  reside  in  triangularized  form  as  individual  columns  in  the  matrix  KELM. 
Dividing  up  the  elements  results  in  a  separate  column  in  KELM  for  in-plane  and  bending 
properties  for  each  superimposed  composite  element.  This  permits  multiplication  of  the 
column  representing  membrane  properties  by  a  linear  element  change  factor  while  the 
column  containing  the  bending  properties  can  be  altered  by  a  nonlinear  change  factor. 
These  manipulations  may  be  done  by  selective  use  of  the  MATMOD  DMAP  module. 

In  the  MSC/NASTRAN  static  solution  procedure,  SOL  24.  each  element  stiffness 
matrix  and  element  mass  matrix  has  its  upper  triangular  and  diagonal  components  stored 
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as  a  column  in  the  assembled  system  matrices  KELM  and  MELM,  respectively.  To 
accomplish  optimal  redesign  it  is  first  necessary  to  checkpoint  these  matrices  and  write 
them  to  FORTRAN  readable  on  disk  using  OUTPUT4.  CBAR  elements  have  element 
matrices  are  already  in  the  basic  coordinate  system;  therefore,  no  coordinate 
transformation  is  necessary  in  order  to  redesign.  It  is  necessary,  however,  to  regenerate 
the  full  symmetric  stiffness  and  mass  matrices  for  each  element.  This  is  done  by  a 
FORTRAN  postprocessor. 

After  the  initial  structural  problem  is  solved,  the  constants  in  the  equation  for  the 

frequency  constraint  must  be  obtained.  First,  the  element  matrices  which  are  output  by 

MSC/NASTRAN  in  diagonalized  form  are  repacked  into  full  matrices  by  FORTRAN 

programming.  Then,  these  matrices  are  used  by  a  DMAP  sequence  to  obtain  the  equation 

components.  This  set  of  DMAP  steps  is  shown  below  for  the  case  with  no  Coriolis  forces. 

ID  GANS,BRONZ3 
TIME  30 
DIAG  4,14 
BEGIN  $ 

FILE  EMAT  =  APPEND  $ 

$  THIS  DMAP  COMPUTES  THE  COMPONENTS  FOR  THE  PREDICTOR  STEP 
$  FOR  THE  BRONZ  PROBLEM. 

$  IT  INVOLVES  8  CBAR  ELEMENTS 

PARAMR  //MPY/V,N,OMEGA2/V,Y,OMEGA/V,Y, OMEGA;  $  SQUARE  OMEGA 
PARAMR  //C,N,COMPLEX//V,N,OMEGA2/O.0/V,N,LAMA/  $  CHANGE 
$  OMEGA2  TO  COMPLEX 
PRTPARM  //  $ 

$  THE  NEXT  SECTION  INVOLVES  THE  CQUAD  ELEMENTS 
LABEL  L2  $ 

INPUTT 4  /K1,M1,D1,UGV l,/4/8/0  $  READ  IN  MATRICES 
ADD  Kl,/K1M/1.0/  $  CUBIC  APPROXIMATIONS  ALREADY  ACCOMPLISHED 
ADD  K1M,D1/K1N/  $  ADD  IN  DIFF  STIFFNESS  FOR  EXTENSIONAL  ELMTS 
SMPYAD  UGVl.KlN.UGV 1. „/A2/3/ 1/1/0/ 1/  $  GET  COMPONENTS  FOR  REDES 
SMPYAD  UGV 1, Ml, UGV1,„/YY 2/3/1/ 1/0/1/  $  GET  COMPONENTS  FOR  REDESIGN 
ADD  YY 2./B2/LAM A'  $  MULTIPLY  YYl  BY  LAMA 

ADD  A2.B2/C2/1.0 /- 1.0/  $  SUBTRACT  MASS  FROM  STIFFNESS  COMPONENTS 
ADD  C2,/E2/C,Y.MINV=  1.0/  $  DIVIDE  BY  GENERALIZED  MASS 
$MATPRN  E2//  $  PRINT  OUT  RESULTS 
APPEND  E2./EMAT/2  $ 

REPT  L2,7  $ 

OUTPUT 4  EMAT ,„,//-  1/8  $ 

MATPRN  EMAT//  $ 

END  $ 

CEND 

TITLE  =  ASSEMBLY  OF  MATRICES  INTO  HOFF’S  EQ.  IV.  19 


SUBTITLE  =  FOR  BRONZ  PROBLEM 
BEGIN  BULK 

PARAM*  OMEGA  2.95200 1E3 

PARAM*  MINV  8.955021E3 

ENDDATA 

When  the  Coriolis  case  is  used,  there  are  two  equations  of  constraint.  One  set 

relates  the  real  components  to  each  other.  These  are  the  magnitude-like  factors.  The 

other  equation  relates  the  imaginary  components  to  each  other.  These  are  the  phase-like 

factors.  The  first  set  of  DMAP  creates  the  real  constraint.  The  second  set  generates  the 

imaginary  constraint. 

ID  GANS,BRONZ45A 
TIME  30 
DIAG  4,14 
BEGIN  $ 

FILE  EMAT  =  APPEND  $ 

$  THIS  DMAP  COMPUTES  THE  COMPONENTS  FOR  THE  PREDICTOR  STEP 
$  FOR  THE  BRONZ  PROBLEM  1ST  EQUATION 
$  IT  INVOLVES  8  CBAR  ELEMENTS 

PARAMR  //MPY/V.N,OMEGA2/V,Y,OMEGA/V,Y, OMEGA/  $  SQUARE  OMEGA 
PARAMR  //C,N,COMPLEX//V,N,OMEGA2/O.0/V,N.LAMA/  $  CHANGE  OMEGA2  TO 
$COMPLEX 

PARAMR  //C,N.COMPLEX//V,Y,OMEGA/O.0/V,N.OMEGAC/  $CHANGE  OMEGA  TO 

$CMPLX 

PRTPARM  //  $ 

$  THE  NEXT  SECTION  INVOLVES  THE  CBAR  ELEMENTS 
LABEL  L2  $ 

INPUTT 4  /Kl,Ml,Dl„/3/8/0  $  READ  IN  MATRICES 
INPUTT 4  /UGVRl,UGVTl,Bl„/3/8/0  $  READ  IN  MATRICES 
MATPRN  Kl,Ml,Dl//  $ 

MATPRN  UGVR1,UGVI1,B1//  $ 

ADD  K1,/K1M/1.0/  $  CUBIC  APPROXIMATIONS  ALREADY  ACCOMPLISHED 
ADD  KlM.Dl/KlN/  $  ADD  IN  DIFF  STIFFNESS  FOR  EXTENSION AL  ELMTS 
SMPYAD  UGVRl,KlN.UGVRl,„/A2/3/l/ 1/0/1/  $  GET  COMPONENTS  FOR  REDES 
SMPYAD  UGVRl.Ml.UGVRl.„/YY2/3/l/l/0/l/  $  GET  COMPONENTS  FOR  REDESIGN 
SMPYAD  UGVIl.Ml.UGVIl,, ./XX2/3/ 1/ 1/0/ 1.  $ 

SMPYAD  UGVIl.KlN.UGVH.„/rWW2/3, 1/1/0/1/  $  GET  COMPONENTS  FOR 
SREDESIGN 

ADD  YY 2  .XX2/Y 21 1.0/—  1.0  $ 

ADD  Y  2, /B  2/LAM  A/  $  MULTIPLY  YYl  BY  LAMA 

ADD  A2.WW2/CC2/1.0/-  1.0/  $  SUBTRACT  MASS  FROM  STIFFNESS  COMPONENTS 
ADD  CC2,B2/E2/1.0/-  1.0/  $  SUBTRACT  CORIOLIS 
APPEND  E2,/EMAT/2  $ 

REPT  L2,7  $ 

INPUTT 4  /UGVR,UGVI,KGGX„/3/8/0  $  READ  IN  DISPL  AND  CORIOLIS  MATRICES 
INPUTT 4  /MGGX,KDNN,BCOR„/3/8/0  $  READ  IN  MORE  MATRICES 
MATPRN  MGGX//  $ 

SMPYAD  UGVI,BCOR,UGVI„,/Tl/3/l/l/0/l/  $  GET  COMPONENTS  FOR  RHS 
SMPYAD  UGVR,MGGX,UGVR,„/VV  1/3/1/ 1/0/1/  $  GET  COMPONENST 


SMPYAD  UGVI,MGGX,UGVI,„/VV 2/3/1/ 1/0/1/  $ 

ADD  W1.W2/V  1/1.0/- 1.0  $ 

MATPRN  Tl.Vl//  $ 

OUTPUT4  EMAT,,,,//- 1/8  $ 

MATPRN  EMAT//  $ 

END  $ 

CEND 

TITLE  =  ASSEMBLY  OF  MATRICES  INTO  GANS  EQUATION  1 
SUBTITLE  =  FOR  BRONZ  PROBLEM  WITH  CORIOLIS,  PREDICTOR 
BEGIN  BULK 

PARAM*  OMEGA  2.951942E3 

ENDDATA 

ID  GANS,BRONZ45B 
TIME  30 
DIAG  4,14 
BEGIN  $ 

FILE  EMAT  =  APPEND  $ 

$  THIS  DMAP  COMPUTES  THE  COMPONENTS  FOR  THE  PREDICTOR  STEP 
$  FOR  THE  BRONZ  PROBLEM  2ND  EQUATION 
$  IT  INVOLVES  8  CBAR  ELEMENTS 

PARAMR  //MPY/V,N,OMEGA2/V,Y.OMEGA/V,Y, OMEGA'  $  SQUARE  OMEGA 
PARAMR  //C,N,COMPLEX//V.N,OMEGA2/O.0/V,N.LAMA/  $  CHANGE  OMEGA2  TO 
ICOMPLEX 

PARAMR  //C,N,COMPLEX//V.Y, OMEGA/0. 0/V,N, OMEGA C/  $CHANGE  OMEGA  TO 

$CMPLX 

PRTPARM  //  $ 

$  THE  NEXT  SECTION  INVOLVES  THE  CBAR  ELEMENTS 
LABEL  L2  $ 

INPUTT 4  /Kl,Ml,Dl„/3/8/0  $  READ  IN  MATRICES 
INPUTT 4  AJGVRl,UGVIl,Bl,,/3/8/0  $  READ  IN  MATRICES 
ADD  Kl,/KlM/1.0/  $  CUBIC  APPROXIMATIONS  ALREADY  ACCOMPLISHED 
ADD  K1M,D1/K1N/  $  ADD  IN  DIFF  STIFFNESS  FOR  EXTENSIONAL  ELMTS 
SMPYAD  UGVRl, KIN, UGVIl,„/A2T/3/ 1/1/0/ 1/  $  GET  COMPONENTS  FOR  REDES 
ADD  A2T,/A2/2.0/  $  DOUBLE  A2T 

SMPYAD  UGVRl, Ml, UGVIl,„/YY2T/3/l/l/0/l/  $  GET  COMPONENTS  FOR 

$REDESIGN 

ADD  YY2T./YY2/2.0/  $ 

SMPYAD  UGVRl, Bl, UGVRl,,, /ZZ2/3/1/1/0/1/  $  GET  COMPONENTS  FOR  REDESIGN 
SMPYAD  UGVIl,Bl.UGVTl„,/XX2/3/l/l/0/l/  $ 

ADD  ZZ2,XX2/ZZ3 

ADD  YY2./B2/LAMA.'  $  MULTIPLY  YYl  BY  LAMA 
ADD  ZZ3,/Z2/OMEGAC/  $  MULTIPLY  ZZ3  BY  OMEGAC 

ADD  A2,B2/CC2/1.0/—  1.0/  $  SUBTRACT  MASS  FROM  STIFFNESS  COMPONENTS 
ADD  CC2.Z2/E2/1. 0/1.0/  $  ADD  CORIOLIS 
SMATPRN  A2.B2.Z2//  $  PRINT  OUT  RESULTS 
APPEND  E2./EMAT/2  $ 

REPT  L2.7  $ 

INPUTT 4  /UGVR.UGVI,KGGX„/3/8/0  $  READ  IN  DISPL  AND  CORIOLIS  MATRICES 
INPUTT 4  /MGGX,KDNN,BCOR,./3/8/0  $  READ  IN  MORE  MATRICES 
SMPYAD  UGVR,BCOR,UGVR„./TT  1/3/1/ 1/0/1/  $  GET  COMPONENTS  FOR  RHS 
SMPYAD  UGVI,BCOR,UGVI,„/TT2/3/l/ 1/0/1/  $ 

ADD  TT1,TT2/Tl  $ 

SMPYAD  UGVR,MGGX,UGVI,„/Wl/3/l/l/0/l/  $ 
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ADD  W1./V1/2.0  $ 

MATPRN  T1,V1//  $  OUTPUT  RESULT 
OUTPUT4  EMAT,,,,//- 1/8  $ 

MATPRN  EMAT//  $ 

END  $ 

CEND 

TITLE = ASSEMBLY  OF  MATRICES  INTO  GANS  EQUATION  2 
SUBTITLE  =  FOR  BRONZ  PROBLEM  WITH  CORIOLIS,  PREDICTOR 
BEGIN  BULK 

PARAM*  OMEGA  2.951942E3 

ENDDATA 

Once  the  equations  of  constraint  are  obtained,  the  problem  in  optimization  is  ready 
to  be  run.  This  is  done  by  the  use  of  a  FORTRAN  program  that  calls  the  ADS  program, 
which  is  the  optimizer  (Vanderplaats  1985).  The  objective  function  can  be  either  minimum 
weight  or  minimum  change.  This  produces  the  values  for  the  element  change  factors  ae 
and  hence  the  thicknesses  of  the  predictor  system.  An  intermediate  finite  element  analysis 
is  then  done  to  generate  new  system  matrices  and  eigenvectors.  This  is  the  conclusion  of 
the  predictor  step.  As  stated  in  Chapter  6,  the  predictor  results  represent  a  linear  solution 
of  the  problem  in  optimization. 

In  the  corrector,  the  MSC/NASTRAN  finite  element  analysis  represented  by  the 

last  step  of  the  predictor  becomes  the  first  step  of  the  corrector.  The  structural  matrices  so 

created  are  used  in  a  corrector  DMAP  to  produce  the  energy  balance  constraint  equation. 

For  the  case  involving  no  Coriolis  forces,  there  is  one  constraint  equation.  When  Coriolis 

effects  are  included,  there  are  two  equations.  The  following  DMAP  represent  the  steps 

necessary  for  the  real  case. 

ID  GANS.BRONZ7 
TIME  30 
DIAG  4,14 
BEGIN  $ 

FILE  EMAT  =  APPEND  $ 

$  THIS  DMAP  COMPUTES  THE  COMPONENTS  FOR  THE  CORRECTOR  STEP 
$  FOR  THE  BRONZ  PROBLEM.  IT  INVOLVES 
$  8  CBAR  ELEMENTS. 

PARAMR  //MPY/V,N,OMEGA2/V,Y.OMEGA/V.Y. OMEGA/  $  SQUARE  OMEGA 
PARAMR  //C,N,COMPLEX//V,N,OMEGA2/0.0/V,N.LAMA/  $  CHANGE  OMEGA2| 

$  TO  COMPLEX 
PRTPARM  //  $ 

$  THE  NEXT  SECTION  INVOLVES  THE  CQUAD  ELEMENTS 
LABEL  L2  $ 
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INPUTT 4  /Kl,Ml,Dl,UGVl,/4/8/0  $  READ  IN  MATRICES 
$MATPRN  K1,M1,D1,UGV1//  $ 

ADD  K1./K1M/1.0/  $  CUBIC  APPROXIMATION  ALREADY  ACCOMPLISHED 
ADD  K1M.D1/K1N/  $  ADD  IN  DIFF  STIFFNESS  FOR  EXTENSIONAL  ELMTS 
SMPYAD  UGV  1,K1N,UGV l„,/A2/3/l/ 1/0/1/  $  GET  COMPONENTS  FOR  REDES 
SMPYAD  UGV  1, Ml, UGV l,„/YY2/3/l/l/0/l/  $  GET  COMPONENTS  FOR  REDESIGN 
ADD  YY2/B2/LAMA/  $  MULTIPLY  YYl  BY  LAMA 

ADD  A2,B2/E2/1.0/— 1.0/  $  SUBTRACT  MASS  FROM  STIFFNESS  COMPONENTS 
$  MATPRN  E2//  $  PRINT  OUT  RESULTS 
APPEND  E2,/EMAT/2  $ 

REPT  L2,7  $ 

$  THE  NEXT  SECTION  COMPUTES  THE  RHS  OF  THE 
$  CORRECTOR  EQUATION 

$  THE  FOLLOWING  MATRICES  ARE  OUTPUT  FROM  THE  WHITE 6 
$  PROBLEM  AND  STUCK  ON  THE  END  OF  THE  RESULTS  OF  REPAK.FOR 
INPUTT4  /UGV, KGGX, MGGX.KDNN.,  4/8/0  $ 

$MATPRN  UGV,KGGX,MGGX,KDNN//  $ 

MATMOD  UGV„„./UGV l,/l/C,Y,COLNUM  =  1  $ 

$MATPRN  UGV1  //  $ 

SMPYAD  UGVl,MGGX,UGVl,„/E3/3////l  $ 

ADD  KGGX,KDNN/KNEW/1. 0/1.0  $ 

EQUTV  KNEW, KGGX' ALWAYS  $ 

$MATPRN  KGGX//  $ 

SMPYAD  UGVl.KGGX,UGVl,„/F/3////l  $ 

ADD  E3, /E/LAMA  $ 

ADD  E,F/G/1.0/~  1.0  $ 

OUTPUT 4  EMAT,,,,//-  1/8  $ 

MATPRN  E,F,G,EMAT//  $ 

END  $ 

CEND 

TITLE  =  ASSEMBLY  OF  MATRICES  INTO  CORRECTOR 
SUBTITLE = FOR  GRAY  PROBLEM.  2ND  10%  CHANGE 
BEGIN  BULK 

PARAM*  OMEGA  3.247201E3 

ENDDATA 

This  next  DMAP  is  the  program  for  running  the  problem  involving  Coriolis  forces. 

ID  GANS,BRONZ49A 
TIME  30 
DIAG  4,14 
BEGIN  $ 

FILE  EMAT  =  APPEND  $ 

$  THIS  DMAP  COMPUTES  THE  COMPONENTS  FOR  THE  CORRECTOR  STEP 
$  FOR  THE  BRONZ  PROBLEM.  IT  INVOLVES 
$  8  CBAR  ELEMENTS. 

PARAMR  //MPY/V.N,OMEGA2/V,Y,OMEGA/V.Y, OMEGA/  $  SQUARE  OMEGA 
PARAMR  //C,N,COMPLEX//V,N,OMEGA2/O.0/V,N,LAMA/  $  CHANGE  OMEGA2  TO 
$COMPLEX 

PARAMR  //C,N,COMPLEX//V,Y,OMEGA/0.0/V,N,OMEGAC/  $  CHANGE  OMEGA 

PRTPARM  //  $ 

$  THE  NEXT  SECTION  INVOLVES  THE  CQUAD  ELEMENTS 
LABEL  L2  $ 

INPUTT 4  /Kl,Ml,Dl„/3/8/0  $  READ  IN  MATRICES 
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INPUTT 4  /UGVRl, UGVIl, Bl„/3/8/0  $  READ  IN  MATRICES 
$MATPRN  K1.M1.D1//  $ 

$MATPRN  UGVRl, UGVIl, Bl//  $ 

ADD  Kl,/K1M/1.0/  $  CUBIC  APPROXIMATION  ALREADY  ACCOMPLISHED 
ADD  KlM,Dl/KlN/  $  ADD  IN  DIFF  STIFFNESS  FOR  EXTENSIONAL  ELMTS 
SMPYAD  UGVRl, KIN, UGVRl,„/A2A/3/l/l/0/l/  $  GET  COMPONENTS  FOR  REDES 
SMPYAD  UGVIl, KIN, UGVIl,„/A2B/3/l/ 1/0/1/  $ 

ADD  A2A,A2B/A2/1.0/-  1.0  $ 

SMPYAD  UGVRl, Ml, UGVRl,,, /YY2A/3/1/ 1/0/1/  $  GET  COMPONENTS  FOR 
$REDESIGN 

SMPYAD  UGVIl, Ml.UGVn,„/YY2B/3/l/l/0/l/  $ 

ADD  YY2A,YY2B/YY2/1.0/- 1.0  $ 

ADD  YY 2 ,/B 2/LAMA/  $  MULTIPLY  YYl  BY  LAMA 

ADD  A2,B2/E2/1.0/— 1.0/  $  SUBTRACT  MASS  FROM  STIFFNESS  COMPONENTS 
$MATPRN  A2,B2,Z2//  $  PRINT  OUT  RESULTS 
APPEND  E2,/EMAT/2  $ 

REPT  L2,7  $ 

$  THE  NEXT  SECTION  COMPUTES  THE  RHS  OF  THE 
$  CORRECTOR  EQUATION 

$  THE  FOLLOWING  MATRICES  ARE  OUTPUT  FROM  THE  WHITE 6 
$  PROBLEM  AND  STUCK  ON  THE  END  OF  THE  RESULTS  OF  REPAK.FOR 
INPUTT4  /UGVR,UGVT,BCOR„/3/8/0  $ 

INPUTT4  /KGGX,MGGX,KDNN„/3/8/0  $ 

3MATPRN  KGGX,MGGX,KDNN,BCOR//  $ 

$MATPRN  UGV1//  $ 

SMPYAD  UGVR,MGGX,UGVR,„/E3A/3/l/l/0/l  $ 

SMPYAD  UGVI,MGGX,UGVI,„/E3B/3/l/l/0/l  $ 

ADD  E3A,E3B/E3/1.0/— 1.0  $ 

ADD  KGGX,KDNN/KNE W / 1 . 0/ 1 . 0  $ 

EQUIV  KNEW, KGGX/ ALWAYS  $ 

SMATPRN  KGGX//  $ 

SMPYAD  UGVR,KGGX,UGVR,„/FA/3/l/l/0/l  $ 

SMPYAD  UGVI,KGGX,UGVI,„/FB/3/ 1/1/0/ 1  $ 

ADD  FA,FB/F/1.0/-  1.0  $ 

ADD  E3, /E/LAMA  $ 

ADD  E,F/I/1.0/—  1.0  $ 

OUTPUT4  EMAT,,,,//-  1/8  $ 

MATPRN  E,F,I,EMAT//  $ 

END  $ 

CEND 

TITLE  =  ASSEMBLY  OF  MATRICES  INTO  CORRECTOR  EQUATION  1 
SUBTITLE  =  FOR  BRONZ  PROBLEM,  MIN  STRUCT  CHANGE,  CORIOLIS 
BEGIN  BULK 

PARAM*  OMEGA  3.247136E3 

ENDDATA 

ID  GANS,BRONZ49B 
TIME  30 
DIAG  4,14 
BEGIN  $ 

FILE  FMAT  =  APPEND  $ 

$  THIS  DMAP  COMPUTES  THE  COMPONENTS  FOR  THE  CORRECTOR  STEP 
$  FOR  THE  BRONZ  PROBLEM.  IT  INVOLVES 
$  8  CBAR  ELEMENTS. 


PARAMR  //MPY/V,N,0MEGA2/V,Y,0MEGA/V,Y, OMEGA/  $  SQUARE  OMEGA 
PARAMR //C,N,COMPLEX//V,N,OMEGA2/0.0/V,N, LAMA/  $  CHANGE  OMEGA2  TO 
$COMPLEX 

PARAMR  //C, N,COMPLEX//V,Y, OMEGA/0. 0/V,N,OMEGAC/  $  CHANGE  OMEGA 
PRTPARM  //  $ 

$  THE  NEXT  SECTION  INVOLVES  THE  CQUAD  ELEMENTS 
LABEL  L2  $ 

INPUTT 4  /Kl,Ml,Dl„/3/8/0  $  READ  IN  MATRICES 
INPUTT4  /UGVRl,UGVIl,Bl„/3/8/0  $  READ  IN  MATRICES 
$MATPRN  K1.M1.D1//  $ 

$MATPRN  UGVR1,UGVI1,B1//  $ 

ADD  K1./K1M/1.0/  $  CUBIC  APPROXIMATION  ALREADY  ACCOMPLISHED 
ADD  K1M.D1/K1N/  $  ADD  IN  DIFF  STIFFNESS  FOR  EXTENSIONAL  ELMTS 
SMPYAD  UGVRl,KlN,UGVTl,„/A2T/3/l/ 1/0/1/  $  GET  COMPONENTS  FOR  REDES 
ADD  A2T./A2/2.0  $ 

SMPYAD  U G VR 1  ,M  1  ,UG VI 1 , , ,/YY 2T/3/ 1/1/0/ 1/  $  GET  COMPONENTS  FOR 

$REDESIGN 

ADD  YY2T./YY 2/2.0  $ 

SMPYAD  UGVRl,Bl,UGVRl,„/ZZ2A/3/l/l/0/l/  $  GET  COMPONENTS  FOR  REDES 
SMPYAD  UGVIl.BlfUGVIl,„/ZZ2B/3/l/ 1/0/1/  $ 

ADD  ZZ2A.ZZ2B/ZZ2/  $ 

ADD  YY2./B2/LAMA J  $  MULTIPLY  YY1  BY  LAMA 
ADD  ZZ2./Z2/OMEGAC/  $  MULTIPLY  ZZ2  BY  OMEGAC 

ADD  A2,B2/C2/1.0/—  1.0/  $  SUBTRACT  MASS  FROM  STIFFNESS  COMPONENTS 
ADD  C2.Z2/E2/1. 0/1.0/  $  ADD  CORIOLIS  TERMS 
$MATPRN  A2,B2,C2,Z2 //  $  PRINT  OUT  RESULTS 
APPEND  E2./FMAT/2  $ 

REPT  L2.7  $ 

$  THE  NEXT  SECTION  COMPUTES  THE  RHS  OF  THE 
$  CORRECTOR  EQUATION 

$  THE  FOLLOWING  MATRICES  ARE  OUTPUT  FROM  THE  WHITE6 
$  PROBLEM  AND  STUCK  ON  THE  END  OF  THE  RESULTS  OF  REPAK.FOR 
INPUTT4  /UGVR,UGVT,BCOR„/3/8/0  $ 

INPUTT 4  /KGGX,MGGX,KDNN„/3/8/0  $ 

MATPRN  KGGX,MGGX,KDNN.BCOR//  $ 

$MATPRN  UGVl//  $ 

SMPYAD  UGVR,MGGX,UGVI,„/E3T/3/l/l/0/l  $ 

ADD  E3T./E3/2.0  $ 

ADD  KGGX.KDNN/KNEW/1.0/1.0  $ 

EQUIV  KNEW, KGGX/ ALWAYS  $ 

$MATPRN  KGGX//  $ 

SMPYAD  UGVR.KGGX.UG VI,.. /FT/3/ 1/ 1/0/1  $ 

ADD  FT./F/2.0  $ 

ADD  E3, /E/LAMA  $ 

SMPYAD  UGVR,BCOR,UGVR,„/H3A/3/l/ 1/0/1  $ 

SMPYAD  UGVl, BCOR,UGVI,„/H3B/3/l/ 1/0/1  $ 

ADD  H3A.H3B/H3/  $ 

ADD  H3./H/OMEGAC  $ 

ADD  E.F/I/1.0/- 1.0  $ 

ADD  I.H/G/1.0/-  1.0  $ 

OUTPUT 4  FMAT,,,,//-  1/8  $ 

MATPRN  E,F,H,G,FMAT//  $ 

END  $ 

CEND 
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TITLE  =  ASSEMBLY  OF  MATRICES  INTO  CORRECTOR  EQUATION  1 
SUBTITLE  =  FOR  BRONZ  PROBLEM,  MIN  STRUCT  CHANGE,  CORIOLIS 
BEGIN  BULK 

PARAM*  OMEGA  3.247136E3 

ENDDATA 

Next,  the  optimization  step  for  the  corrector  is  run.  The  objective  function  in  this 
part  is  either  minimum  change  or  minimum  weight. 

Once  the  corrected  values  of  ae  are  obtained,  the  element  thicknesses  for  the  system 
can  be  calculated.  Then,  another  finite  element  analysis  is  run.  This  is  required  for  the 
intermediate  iterative  or  incremental  steps,  but  otherwise  it  is  only  used  to  check  the  final, 
optimized  solution. 


APPENDIX  B 


BEAM  SOLUTIONS 

Derivation  of  Exact  Solution  for  Beam  Vibration 

In  this  Appendix,  the  exact  solution  for  the  nonrotating  beam  problem  posed  in 
Chapter  4  will  be  obtained.  Figure  29  illustrates  the  problem  being  considered:  it  is  the 
planar  problem  of  a  (nonrotating)  beam  of  constant  section  with  one  end  clamped  and  the 
other  end  mounted  on  a  roller  that  permits  only  vertical  (y-axis)  displacements  and  no 
rotations.  The  beam  has  modulus  of  elasticity  E,  moment  of  inertia  I,  length  L,  linear 
mass  density  m,  density  p,  and  cross-sectional  area  A.  Let  u  be  the  eigenfrequency  of  the 
structure  and  Y(x)  denote  the  vertical  displacement  as  a  function  of  beam  wise  distance  x. 


y 

I 


Figure  29.  Rotating  Beam 


The  governing  differential  equation  for  this  problem  is: 


d2  [  d2Y(x)  1 

-  El -  =w2mY(x) 

dx2  dx2 


The  boundary  conditions  are: 


Y(0)  =  0 
Y'(0)  =  0 
Y'(L)  =  0 
Y”(L)  =  0 


Defining: 


g*  =  i£E. 

“  El 

and  noting  that  both  E  and  I  are  constant  in  this  problem  obtains  a  new  version  of  the 


differential  equation: 

d4Y(x) 


-  /T4 Y(x)  =  0 


The  solution  to  the  above  equation  may  be  assumed  to  be  of  the  form: 

Y  (x)  =  C  jsin/Jx  +  C2cosi3x  +  C3sinh/3x  +  C4coshSx  (B.  4) 

Taking  derivatives  of  Equation  (B.4)  and  applying  the  boundary  conditions  results  in 
the  following  matrix  equation  where  the  unknowns  are  the  Cj  constants. 


1  u  1  u  u2  -  ;rn 

cosSL  —  sin/3L  ccsh<3L  sinh5L  C3 

-cosBL  sin5L  cosh<3L  sinh^L  C4 

For  Equation  (B.5)  to  have  a  nontrivial  solution,  the  matrix  on  the  left-hand  side  of  the 
equation  must  have  zero  determinant.  Solving  for  the  first  mode  one  obtains: 

0L  =  2.4259  0 

Substituting  for  3,  one  obtains  the  following  expression  for  w: 
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u/=  5.8850 


El 


mL4 


(1/2) 


(B.7) 


For  this  problem,  the  beam  is  made  of  IN-718  steel,  with  E  of  73.77E3  MPa,  I  equal 
to  3.2252E4  mm4,  A  equal  to  625  mm2,  p  of  2.774E-9  Mg/mm3,  L  equal  to  250  mm,  and 
m  calculated  to  be  1.73375E-6  Mg/mm.  Using  these  values,  the  fundamental 
eigenfrequency  is  given  as: 

Wj  =  3.5043E3  rad/sec 
=  (theoretical) 

This  is  the  value  used  to  calibrate  the  finite  element  model  for  the  rotating  beam. 
Convergence  of  as  determined  by  MSC/NASTRAN  finite  element  analysis  to  the 
theoretical  value,  as  the  model  is  refined,  is  shown  in  Chapter  6,  Table  1. 

Derivation  of  Tangent  Stiffness  Matrix 

The  nonlinear  problem  under  consideration  in  this  dissertation  is  called  geometric 
nonlinearity  and  the  equations  of  equilibrium  must  be  reformulated  for  the  deformed 
configuration.  An  incremental  procedure  can  then  be  used  to  obtain  a  tangent  stiffness 
matrix  (Przemieniecki,  1985). 

For  a  linear  elastic  material,  the  relationship  between  stress  {a-}  and  strain  {e}  is 
given  by: 

{cr}  =  [G]{e}  (B.8) 

where  [G]  is  the  matrix  that  relates  strain  to  stress.  Since  the  material  is  linear,  the  same 
matrix  [G]  relates  the  change  in  strain  {Ac}  to  the  change  in  stress  {Ac}: 

{Act}  =  [G]{Ac}  (B.91 

Since  the  problem  is  geometrically  nonlinear,  one  matrix  [Bs]  is  needed  to  relate 
displacement  {u}  to  strain  {«}  and  another  matrix  [Bl]  is  required  to  related  the  change  in 
displacement  {Au}  to  the  change  in  strain  {Ae}  (Anderson,  1985): 

(B.  10) 


1 


V 


«  =  CBs]{u} 
{Ac}  =  [Bl]{Au} 


(B.ll) 


3rTC  V  VXTWWH7>.7  KT. 


!  Kwmw VTWVJVm  vy  VT 


n'^vr^^'f.^rew  ^Tvrr»'4»v"  ..^ 


107 


The  matrix  [Bl]  is  at  most  linear  and  can  be  taken  to  be  the  sum  of  the  matrix  that 
relates  displacements  to  strains  for  infinitesimally  small  displacements  [Bq]  and  a 
transformation  matrix  linearly  dependent  on  displacements: 

[Bl]  =  [Bq]  +  [B\(u)]  (B.  12) 

The  tangent  stiffness  matrix  [Klj  can  be  found  by  incrementing  the  applied  load  f 
and  the  displacement  u  such  that  the  equation  of  equilibrium  for  the  structure  becomes: 

A(/v[Bl]TMdV)  =  A{f}  (B.  13) 

Expanding  the  left-hand  side  of  the  above  equation  one  obtains: 

/v[ABl]T{<r}dV  +  /v  [Bl]T{Acr}dV  =  {Af}  (B.14) 

The  first  integral  gives  the  differential  stiffness  matrix: 

[KD]{Au}  =  /v[ABl]V}dV  (B.  15) 


An  alternative  definition  for  the  differential  stiffness  matrix  is  given  by  Cook  (1974)  and  is 
shown  in  Equation  (5.1). 

Equations  (B.9),  (B.ll),  and  (B.  12)  are  applied  to  the  second  integral  of  Equation 
(B.14)  so  that  Equation  (B.14)  may  be  expressed  as: 

[Kx]{Au}  =  {Af}  (B.  16) 

where 


[K0] 

fKT]  =  [KD]  +  /v[B‘]TfG][B^dV 

+  fv  [Bc>]T[G][Bj]dV  +  /v  rBt1]T[G][B'l]dV  +  /v  [Bt1]T[G][Bt1]dV 

- - -  .  _  v 

[Kl] 


(B.  17"! 


The  first  integral  on  the  right  hand  side  of  Equation  (B.17)  is  the  conventional, 
infinitesimally  small  displacement  stiffness  matrix,  [K0J.  The  other  three  integrals 
together  comprise  the  1"_  e  displacement  stiffness  matrix  [KL],  Therefore,  Equation 
(B.17)  can  be  rewritten  as: 


1 


a 


. 

\ 


[Kt]  =  [K0]  +  [Kd]  +  [Kl] 


(B.  18t 


The  above  equation  is  identical  to  Equation  (4.1).  In  this  dissertation,  it  will  be  assumed 
that  [Kl]  can  be  neglected  in  comparison  to  [KoJ  and  [KD], 

As  an  example  of  a  differential  stiffness  matrix,  let  us  examine  the  rotating  bar 
shown  in  Figure  6  in  Chapter  4,  taken  as  a  single  element.  First,  define  three  degrees  of 
freedom  at  each  end:  x-displacement,  z-displacement,  and  y-rotation.  The  applied  stress 
due  to  centrifugal  loading  is  ax  and  is  given  by: 

<7X  =  P/A  (B.  19) 

where  A  is  the  cross-sectional  area  of  the  bar  and  P  is  the  centrifugal  loading  and  is  shown 
in  Equation  (4.19)  to  be  4ir2fl2mL.  Therefore,  the  differential  matrix  can  be  shown  to  be: 


[Kn]  =  f  - 


0  0  0  0  0  0 

0  6/5  L/10  0  -6/5  L/10 

0  L/10  2L2/15  0  -L/10  -L2/30 

0  0  0  0  0  0 

0  -6/5  -L/10  0  6/5  -L/10 

0  L/10  — L2/30  0  -L/10  2L2/15 


(B.20) 
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APPENDIX  C 

DEGRADED  SOLUTIONS  FOR  BLADE  AND  BEAM 

In  this  section,  the  problems  of  the  rotating  blade  and  beam  presented  in  Chapter  7 
will  be  degraded  (simplified  by  omitting  certain  features).  This  will  permit  the  examination 
of  the  results  when  centrifugal  effects  are  neglected  in  the  optimization  and  in  the 
structural  analysis. 

Results  of  Degraded  Cases 

Table  9  shows  the  results  in  optimization  for  the  degraded  problems.  This  table  is 
similar  to  Table  3.  Case  12  is  the  rotating  blade  shown  in  Figure  16  with  no  centrifugal 
effects  in  the  optimization  constraints.  Minimum  change  is  used  both  in  the  predictor  and 
corrector.  Case  13  is  similar  to  Case  12  with  the  hybrid  objective  function  employed:  that 
is,  minimum  weight  is  used  in  the  predictor  and  minimum  change  is  used  in  the  corrector. 
In  Case  14,  centrifugal  effects  are  neglected  both  in  the  optimization  and  in  the  structural 
analysis.  Minimum  change  is  used  in  both  the  predictor  and  the  corrector.  Case  15  also 
neglects  the  centrifugal  forces,  but  the  hybrid  optimization  routine  is  used.  In  Case  16  and 
Case  17,  the  beam  problem  shown  in  Figure  10  is  reanalyzed  with  centrifugal  effects 
negelected  in  the  optimization  but  included  in  the  structural  analysis.  Case  16  uses 
minimum  structural  change  while  Case  17  uses  the  hybrid  optimization. 

Table  10  contains  the  results  of  the  parameters  in  the  optimization  routines.  This 
table  is  similar  to  Table  6. 

Observations 

In  the  blade  optimization,  the  change-change  objective  function  obtains  satisfactory 
results  for  the  frequency  change  even  if  centrifugal  effects  are  completely  ignored. 
However,  there  is  a  slight  amount  of  overshoot  on  the  linear  predictor  step.  The  hybrid 
objective  function  obtains  good  results  at  the  end  of  a  complete  cycle  when  centrifugal 
effects  are  retained  in  the  structural  analysis  (Case  13),  though  the  linear  predictor  does 
not  provide  a  good  solution.  When  centrifugal  effects  are  completely  ignored  (Cases  14  and 
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15),  satisfactory  solutions  are  still  obtained. 

In  the  beam  optimization,  good  results  are  observed  in  Case  16.  However.  Case  17 
exhibits  the  worst  results.  Not  even  the  corrector  can  obtain  the  frequency  goal  without  a 
1  significant  error. 

The  results  shown  in  Table  10  are  similar  to  those  shown  in  the  cases  where 

i 

■  centrifugal  effects  are  included  (see  Tables  6  and  8).  Once  again,  the  optimizer  gives  a 

value  of  the  objective  function  close 
minimum  weight  objective  function. 

’  constraint  is  satisfied. 


I 
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to  zero  for  minimum  change  but  a  negative  number  for 
Values  of  the  constraint  close  to  zero  imply  that  the 
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